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ABSTRACT 


BENNETT, ROBERT MERRILL. An Analytical Investigation of an Oscillating 
Wedge in a Supersonic Perfect Gas Flow. (Under the direction of ROBERT 
W. TRUITT). 

Several aspects of the oscillating wedge are Investigated to evalu- 
ate both the resulting trends for the wedge and methods of analyzing 
unsteady flows. 

The existing hypersonic small disturbance theory for the oscillat- 
ing wedge is considered. Algebraic expressions are obtained for the low 
reduced frequency aerodynamic forces. Series expansions of the result- 
ing forces in powers of l/MG^ reveal weak similarity parameters for the 
pitching rate derivatives that indicate that a simple scaling of flutter 
with ratio of specific heats is not apparent. The importance of the 
reflections of the waves generated by the surface from the shock is 
reemphasized. Aerodynamic damping terms are found to have significant 
influence on flutter even for large values of the mass ratio parameter. 
The corresponding effects of angle of attack are shown to be small for 
the thin wedge. 

The equations are developed for linearized perturbations about the 
known steady flow conditions that can be applied for any wedge angle or 
Mach number for which the local flow is supersonic. The method of 
integral relations is applied to the resulting equations in a one-strip 
approximation. Algebraic expressions for low reduced- frequency aero- 
dynamics are obtained and a set of ordinary differential equations are 
obtained for general oscillatory motion. The method gives accurate 



results for low reduced- frequency. However, for cases in which the 
aerodynamic forces vary rapidly with frequency the results are qualita- 
tively correct but of limited accuracy. Calculations indicate that for 
a range of inclination angles near shock detachment such that the flow 
in the shock layer is transonic , the aerodynamic forces vary rapidly 
both with inclination angle and with reduced frequency indicating 
potentially serious flutter problems near detachment - 

The governing nonlinear flow equations are formulated in a body- 
fixed axis system. The first order numerical finite-difference method 
of Lax is applied and discussed with simplified treatments of the 
boundaries. Although the formulation is suitable for numerical treat- 
ment, and may have merit for unsteady blunt body investigations, the 
need is indicated for further refinement of the treatment of the bound- 
ary conditions and for use of a second order difference scheme. A brief 
quasi-static analysis of oscillating wedge flows indicated that signifi- 
cant nonlinear effects may exist near detachment and for large amplitude 
motions of a thin wedge at high Mach numbers. 
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1. INTRODUCTION 


Aeroelastic problems such as flutter have been a design considera- 
tion for aerospace vehicles for many years, and as such, have been the 
subject of considerable research. The need still exists for further 
refined analytical techniques as the potential economic savings of 
careful design of advanced vehicles can be significant. Recent develop- 
ments in computer technology have led to several highly- developed aero- 
dynamic methods for the analysis and computation of steady flow fields. 
The corresponding unsteady aerodynamics for analyzing flutter, loads, 
and stability have not been as well developed. In this thesis some of 
the techniques that have been applied to the calculation of steady flow 
fields are applied to a simplified configuration — an oscillating two- 
dimensional wedge in a supersonic or hypersonic flow of a perfect gas. 

The steady supersonic flow over a sharp two-dimensional wedge is 
one of the few exact solutions of the nonlinear governing equations of 
inviscid flow. No such corresponding solution for unsteady flow is 
available. However, the uniform steady wedge flow permits significant 
simplification of perturbation methods. This fact has been used to 
advantage in several phenomenological investigations, perturbation 
analyses, and approximating techniques in order to evaluate such 
processes as unsteady wave reflections. Herein the oscillating wedge is 
treated both to gain insight into the application of several methods of 
analysis and computation, and to examine the resulting trends for the 
wedge airfoil. 

First, after discussion of the problem and literature, an existing 
hypersonic small distrubance theory for an oscillating thin wedge is 



2 

extended and applied. A perturbation method involving linearization 
about the known steady flow is then derived and discussed. Subsequently , 
a finite difference technique for calculating the complete unsteady flow 
field of the wedge in motion is presented and discussed in conjunction 
with some calculated quasi-static nonlinear trends. Emphasis herein is 
given to the high supersonic and hypersonic Mach number range and to the 
effects of the ratio of specific heats. 
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2. DISCUSSION OF PROBLEM AND RELATED LITERATURE 

The literature of the general area of aeroelasticity is extensive 
and has been recently reviewed, for example, by Ashley (1970)* Various 
aspects of aeroelasticity have also been discussed by Bisplinghoff and 
Ashley (1962) and Garrick (1969). Supersonic-hypersonic unsteady 
problems have also been reviewed by Wood (1966) and Ashley and Zartarian 
(1961). Thus, the discussion here will be principally concerned with 
several aspects of the specific problem treated in this thesis. 

2.1 The Oscillating Wedge 

The effect of airfoil thickness on flutter at supersonic-to- 
hypersonic speeds can be large and destabilizing (for example, see 
Ashley and Zartarian, 1956; Morgan, et al., 1958; and Runyan and Morgan, 
1962). This effect generally results primarily from a forward chordwise 
shift of the aerodynamic center with increasing thickness. Other 
factors such as lift-curve-slope and pitch damping are also involved, 
however. In order to assess the role of thickness before the destabi- 
lizing effect of thickness was well known, the wedge airfoil was 
attractive for further treatment since the exact steady flow field 
was known. In more recent years, it has also been a basic check case 
for evaluating simplified theories and for phenomenological studies. 
Although the wedge airfoil is often studied for research purposes only, 
such airfoils are infrequently used such as on the vertical tail of the 
X-15 airplane. 

Carrier (1949b) presented a perturbation analysis of the oscillat- 
ing wedge in supersonic flow as based on his earlier study of the 
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stability of the strong and weak shock solutions for steady wedge flow 
(Carrier, 1949a). The analysis was further developed and applied by 
Van Dyke (1953), who gave some brief results. It had been shown (e.g., 
Garrick and Rubinow, 1946) from linearized supersonic flow theory, that 
unstable values of damping- in-pitch of a flat plate exist for forward 
locations of the pitch axis and for Mach numbers less than Vi- 
Van Dyke (1953) demonstrated that the range of Mach numbers and of pitch 
axes for unstable pitch-damping was enlarged for a thin wedge. 

The results of the above-mentioned analysis were used by Van Dyke 
(1954a) to verify the low- frequency results obtained from second order 
theory. Van Dyke (1954a) also presented a solution obtained from second- 
order theory for a thin wedge oscillating at arbitrary frequency. A 
sample comparison with linear theory indicated that the effects of 
thickness on the damping-in-pitch (for k 0) are alleviated at higher 
values of k. 

With the development of the relatively simple piston theory, which 

is applicable to thin surfaces and high Mach numbers, Hayes (1947), 

Lighthill (1953), and Ashley and Zartarian (1956), the wedge airfoil 

was again treated (e.g., Chawla, 1958). Piston theory is limited to 

values of M9 <1 and it has also been noted that stable values of 
w 

the damping-in-pitch are inherent in the piston approximation. 

McIntosh (1965a , 1965b) presented a perturbation solution for the 
oscillating wedge as based on hypersonic small disturbance theory 
(Van Dyke, 1954b), that is essentially a generalization of piston theory 
to include the effects of the shock wave on the local external flow 
properties and the effects of reflections of surface motion-generated 



waves from the shock wave. McIntosh's analysis indicated that inclusion 
of the wave reflection process was essential} otherwise the perturbation 
in pressure approaches infinity as 7 -» 1. Such was observed by Miles 
(i960) in a local flow perturbation analysis. The wave reflection 
process had been earlier studied by Chu (1952). An analysis similar to 
McIntosh's was presented by Appleton (1964) and discussed by Orlik- 
Riickemann (1966b, 1969). The results of McIntosh (1965a) have been 
applied to pitch-plunge flutter of wedge airfoils by Bailie and McFeely 
(1966) and to panels mounted on wedges by Bailie and McFeely (1968). 

The hypersonic small disturbance theory for the oscillating wedge of 
McIntosh (1965a) is again applied and discussed in this thesis in 
Chapter 3* 

Hui (1969a) developed an exact perturbation theory for pitching 
motions of a symmetrical wedge for angles up to detachment and has also 
applied the method to the caret wing. Unstable values of damping- in- 
pitch for the wedge were shown to exist at any Mach number if the wedge 
angle was sufficiently large. For the caret wing, this result was 
eliminated by three-dimensional effects. With an approximate perturba- 
tion solution, Hui (1969b) also indicated that waves generated by the 
moving shock wave can be significant for large wedge angles. A perturba 
tion method is developed here in Chapter 4 that is quite similar to 
that of Hui; however, the method of solution is quite different and the 

O 

analysis is done in more general terms such that any rigid-body motion 
of an inclined- flat compression surface can be described. 

Kuiken (1969) has presented a higher-order perturbation solution 
based on hypersonic small disturbance theory. The quasi-static wedge 
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flow , that is based on the effective wedge angle at the vertex including 
motion, was perturbed to give an effective phase shift of the quasi- 
static solution. The surface pressure was calculated for constant- 
amplitude harmonic motion. As the quasi-static solution is nonlinear, 
the res ul ting pressure waves are unsymmetrical, indicating the presence 
of higher harmonics, and a shift in mean pressure level was noted. 

Curve fits to Kuiken’s results have been given by Orlik-Riickemann 
(1969). The quasi-static solution will be subsequently discussed in 
Chapter 5. In a similar analysis, Kacprzynski (1968) has treated the 
full equations of motion in an approximate solution giving similar 
results. Hui (1970) has subsequently given an exact solution with no 
further numerical results. 

Several other simplified theoretical developments have been given 
that could be applied to the wedge airfoil such as Newtonian flow theory 
(i.e., Zartarian and Sauerwein, 1952). Some have been examined by 
applying them to a diamond airfoil by Yates and Bennett ( 1971) • 

An important aspect of hypersonic flow, at least, is the inter- 
action of the boundary layer with the steady flow field. Some research 
into unsteady effects of a laminar boundary layer on wedges in the weak 
interaction regime is given by King (1966), Orlik-Riickemann (1966a, 
1970), and Hui and East (1971). The general trend is to shift the 
moment center forward and to reduce the damping- in-pitch for forward 
locations of pitch axis. 

Several investigators have also given experimental data relating 
to the overall aerodynamic forces for the pitching wedge. Pugh and 
Woodgate (1963) give the results of two-dimensional tests at supersonic 
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speeds in air. East (1962) has presented some results of measurements 
in air at M = 10, hut with some tip effects. Orlik-Ruckemann (1970) 
has also given some hypersonic results measured in helium. Martuccelli 
(1958) has also presented some measured results at low supersonic speeds. 
Unfortunately , the author is unaware of corresponding two-dimensional 
flutter data as available data are generally for swept-tapered wings 
(i.e., Hanson, 1961). 


2.2 Flow Field Calculations 

Several numerical methods have recently been developed which permit 
computation of highly- comp lex, steady flow fields over various types of 
bodies using the governing nonlinear flow equations and a modem digital 
computer. The proceedings of a recent symposium (NASA publ., 197°) 
gives a broad survey of several of the available techniques. An 
annotated bibliography has also been given by Harlow (1969) and a survey 
volume given by Chu (1968). These methods generally consist of (l) 
direct solution by constructing finite difference approximations to the 
partial differential equations and then solving the resulting system of 
equations numerically, or (2) reducing the governing partial differ- 
ential equations in some approximate manner and then integrating the 
ordinary differential equations. 

Methods of the first type as based on a characteristics approach 
have been given, for example, by Sauerwine (1964, 1965, 1966) and also 
by Rakich (1969). Various other finite-difference approaches have been 
developed and applied, for example, by Barnwell (1971); MacCormack 
(1969), Moretti and Abbett (1966), Lax and Wendroff (196k), and Lax ( 195k). 
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Many of these methods involve determining the steady flow field by the 
time evolution from an assumed approximate initial solution. 

Examples of methods of the second type are the method of integral 
relations (see for example, South, 1968, and Belotserkovskiy and 
Chushkin, 1962) and the method of lines (Klunker, et al., 1971 > for 
example). These methods usually treat the boundary or initial value 
problem directly or iteratively in contrast to the time-asymtotic 
approach . 

Very little work has been done in relation to applying any of these 
methods to unsteady flows except for some simplified one-dimensional 
problems such as given by Bohachevsky and Rubin (19 66), although they 
have been occasionally discussed (i.e., Ashley and Zartarian, 1961) . 

In this thesis two of the above methods are applied to the unsteady 
aerodynamics of the oscillating wedge, the method of integral relations 
is applied to the perturbation equations (Chapter 4) and a first-order 
difference method of Lax (1954) is applied to the full nonlinear 
unsteady flow equations written in a body-axis system (Chapter 5). 

2.5 Comments on Nonlinear Aeroelastie Analysis 

Aeroelastic problems have generally been treated with linearized 
analyses that essentially consider initial tendencies for infinitesimal 
amplitudes. For example, flutter is conventionally posed as a linear 
eigenvalue problem based on linear representations of both struture and 
aerodynamics. The resulting eigenvalues are then independent of 
initial conditions or initial disturbance by the assumption of linearity. 
A no n! i near problem requires considerably more effort to determine the 
overall characteristics as the stability characteristics are dependent 
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on the level of input, etc. Furthermore, the nonlinear characteristics 
are not generally as well known as the level of effort of testing or 
analysis is considerably increased. However, some aeroelastic analyses 
have been performed that consider nonlinearities. Some simplified 
structural nonlinearities have been considered by Woolston, et al. ( 1955 ) 
both experimentally and by analog computation. Nonlinear structural 
representations of panels in the problem of panel flutter have recently 
been considered by Dowell ( 1970 ) and Morino, et al. (1969)* 

The corresponding aerodynamic nonlinear considerations are also 
sparse. However, Landa and Qtrelkov (1963) have given some results of 
an analog computation using Newtonian aerodynamics at high angle of 
attack and including the second harmonic resulting from the derivative 
of C p = 2 sin 2 8. A second dip in the flutter speed was demonstrated 
as o^/o^ -* 2 in addition to the usual dip for -» 1 . The implica- 

tions of such a trend are potentially serious in a multimodal situation 
if such nonlinear effects in the aerodynamics generate large second or 
higher harmonics. Zartarian and Sauerwine (1962) have indicated the 
presence of second harmonics for large amplitude oscillations of a 
biconvex airfoil as calculated from third order piston theory and from 
a local shock-expansion theory, but no application was made. As 
previously discussed the analyses of Kuiken (1968), Kacprzynski (1968), 
and of Hui ( 1970 ) give nonlinear aerodynamic forces on a wedge executing 
simple harmonic motion of low frequency. The quasi-static solution is 
briefly analyzed in Chapter 5 in order to determine where nonlinear 
effects might need consideration for the oscillating wedge. The 
analysis of the quasi-static results is discussed with the results of 
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the finite-difference calculation which has the potential for calculating 
nonlinear forces. It might also be noted that for second order piston 
theory (small M8 W ) , the second harmonic cancels (Lighthill, 1953)* Such 
would not be the case for a flat plate at high angle of attack, for 
example, where only one surface would be effective. 

It might also be noted that Muhlstein ( 1971) has measured large 
second harmonics in the pressure distribution over a wavy wall in the 
transonic range. 



3- A STUDY OF HYPERSONIC SMALL DISTURBANCE THEORY 
FOR AN OSCILLATING WEDGE 
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3.1 Aerodynamic Forces for Pitching and Plunging 

Oscillations 

A thin, two-dimensional wedge undergoing small oscillations in a 
perfect gas flow has been treated by McIntosh (1965a, 1965b) from the 
standpoint of the hypersonic small disturbance theory given by Van Dyke 
(I95 i *-b) . McIntosh's solution differs from piston theory, which is also 
based on hypersonic small disturbance theory, in that it considers the 
change in the local external flow field resulting from the bow shock 
wave; considers the effects of reflection of acoustic waves generated by 
the surface motion from the bow shock onto the surface; and is not limited 
to M9 w « 1. The resulting surface pressure coefficient for linearized 
oscillatory motion about the nonlinear steady flow field (in somewhat 
different notation) is: 


2FK 

C p( X ' t ) = If - 


f;(x) + 2ik f w (x) 


n-1 


2ik(r n -l) 


r (I^x) + 2ik f w (r n x) 


>6 e 


2ikt 


3-1 


where f (x) and f'(x) are the normalized wedge deformation mode shape 


w 


amplitude (f /c), and slope; £ is modal amplitude, h or 0 here; 

k is reduced frequency, ccoy&V^; K = Mt where M is Mach number and 

P , . - 1 -/ 2 

27K; -(7-1) 


t is steady shock wave slope; F = 


2 + (7 - 1 ) K? 


; "K is the 


attenuation factor for disturbances reflected from the shock; and r is 
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a reflection length coefficient- The physical interpretation of (3*1) 

is that the pressure at point x is proportional to the effective local 

slope at x, (f^ + 2ik f^), multiplied by the linearized slope-pressure 
2 

relation — and the factor that gives the local flow correction, FK t ; 
plus the additional portion given by the sum which results from acoustic 
reflection effects, made up of the reflection attenuation (-"\) n , a phase 
factor e 2ik ( r and the effective slope at point (l^x). Thus the 

point x is affected by points upstream Tx, Rx, ..., (as r < l) . 

The surface pressure, coefficient is integrated over the chord of 
airfoil. The resulting aerodynamic forces are expressed as coefficients 
in the form given by Garrick and Rubinow (1946) and subsequently dis- 
cussed and defined herein in section 4-3-3 • The resulting expressions 
for the coefficients for pitching about the leading edge and for plunging 
motions are (multiplied by factors of M and k, and with moments taken 
about the vertex) : 


ML 1 


- Y, <-*) n " r ' 1 ) 

n=l 


(l - cos 0 ) 
' n 



kMLg = FK t 


£ 

n=l 


1+2 ) (-A) 


sin p 
, n 'n 


n 




w 

[ = 8fk t £ (-TO^r 11 - l) 


n sin p 
cos P - n 

n 


0 


n 


n=l 


n 


kMMg 


- 2kML 2 



4 



U) n 


(l - COS P ) 
n 
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k 2 ^ = kML 2 + 8k 2 FK i •I (-A) n (r n )(r n 

n=l 


i) 


cos 3 - S ^ n ^n 

n T - 

- n _ 


n 


kML£ = 


- (MI^) + fk t 


1 + 4 


£ (-A) n (r n ) 

n=l 


cos 3 n -l sin 3 n] 


n 


n 


k 2 MMi = kMMl 
3 2 


l6k 2 FK T £ (-aKi^Ki^-i) 


2 sin 3 2(l - cos 3 ) 

cos 3 S- + n 


n 3 


n 


Pn" 


n=l 


n 


kMM^ = - (MMp 



(3.2) 

The transfer of the moment coefficients to axes other than the leading 
edge is discussed in section 8.J. 

As 3 = 2k (r 11 - l), and A and r are functions of M9 and y 

n W 

only, the coefficients written in the form of equations (3.2) are func- 
tions of 7, M9 , and k, and are valid for M » 1, 9 « 1, 

0 < M9^ < oo, and k < 1 (for the restriction on k see Hayes and 
Probstein (1966, pp. 113-116)). These coefficients have also been given 
by Bailie and McFeely (1966) and McIntosh (1965b) in somewhat different 
form. The form given herein is such that each term should be finite as 
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k ->0. However, note that as k -* 0, 0 n -» 0 and each of expressions 
( 3 * 2 ) involve indeterminant forms and must he given special attention 
for the small values of k of interest in the hypersonic speed range. 
Computational forms for k « 1 are developed by using series expansions 
for the sine and cosine converting the summations of ( 3 * 2 ) to double 
summations. For example ML^ becomes, 

00 “ B 

= 4FK T £(-70 n (l - I* 1 ) £(-l ) m+1 

n=l m=l 

The expressions for K^, k, and r from McIntosh ( 1965 a) are func- 
tions of K = M0 and 7 as follows, 
w 

K = B +Vl + B 2 

T 

where, 

B = l2L^X}JC 


where 


and 


2(K 2 + 1) 

1 _ 1 

(7 + i)k t 2 

4 

D “ (7 + 1 )F 



(3*M 
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vhere ~ l/2 

27K - (7 - 1) 

F = ? 

2 + (7 - 1) K 


Note that the expression for K as given above is that of thin shock 
theory (for example, see Truitt (1959 ) > eqs. 2-155 and 2-156). 


3-2 Low Reduced-Frequency Limiting Case 
In the low frequency (k) expansions for the aerodynamic forces (3*2), 
each inner sum for the sine and cosine expansions is of the form of the 
series of (3«3) involving only even powers of 3 n = 2k (r 11 - l) . If one 
neglects the terms of 0(k ) and higher and retains only the constant 
term, the series over n sum in closed form. The following simplified 
expressions are then obtained. 



(5-5) 


Only the three coefficients ML^, kML^, and kML^ are required 
with (8.12) for a complete description of aerodynamic forces for small 
values of k. As subsequently discussed in Chapter 4, computations have 
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also shown that the effect of frequency is negligible for the applicable 
range of k < 1 except for a k 2 - variation in k 2 ^ and k^Mj for 
7 < 1.2. Thus, the above results can be considered the complete solu- 
tion for the normal range of reduced frequencies for hypersonic flutter. 
Some of the above relations, in different notation, and some approximat- 
ing functions have also been given by Orlik-Ruckemann (1966b and 1969) • 
Substituting for F, and K from equations (3.4) gives 

(3.6) 

2 V(7 + l) 2 K 2 + 16 

which gives a result that can be obtained by differentiating the well 
known thin shock layer result, that is, equation (2-156) of Truitt (1959) . 
Thus, a closed form solution is obtained for the major portion of the 
aerodynamic forces. However, no such reduction has been obtained for 
ML 1 or kML^ . 

It might also be noted that the flutter derivatives used herein, 
and stability derivatives are related for k -» 0 and to first order in 
k as (for example, see Van Dyke (1954a)): 


a 


a 


= 4KL 2 

C m - " 2kM 2 
a 

- -4ii 

0 

P^ 

11 

+ C 7 = 4kL)' 

C + C 

l 4 

q 

m. m 

a q 


(5.7) 


3.3 Limiting Case for Large M9 

w 

For M9 w » 1 for which the effects of wave reflections become 
more significant, limiting forms of (3.5) can be obtained that reveal 
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the role of the shock reflection process and 7 . Expanding F, K t , \ 
and T as series in l/K and substituting in (3-5) gives 


kML_ = (7 + 1) K + b( l/K 5 ) 


(3.8) 


MI^ = 


(r + 1) (2 - 7) K + ( ? + l 1 */- 7? g ) 1, o(l/K 5 ) 

2 WVtT (r + l)(2r _ 1)2 K °W K 1 


ms - K - R - a - t I - o* 1 /* 3 ) 

4 (27 “ 1} (7 + l) (27 - l) 2 K 


+ ™4 ' fay 




(3 .. + . i>y i + 0 (i/k 5 ) 

(7 + l) (27 - 1 r 


The result for kMLg is given by the well-known limiting case of thin 

shock layer theory (Truitt, 1959> eq. 2-157) for large M0 W with 

Cp = (7 + l) 9 w « The parameter (7 + l)K is thus a similarity parameter 

for static pressure and thus for lift-curve slope. The above expressions 

are the consistent relations for the stability derivatives to the same 

level of approximation and the leading terms are identified as similarity 

parameters. Results from the complete expressions (3*5) are compared 

in figure 3.1 with the results from (3-8) retaining terms 0(K) and also 

retaining terms 0(l/K) . The use of the leading terms of (3*8) as 

similarity parameters gives results comparable to the results shown in 

figure 3*1 for retaining only M0 terms. These similarity parameters 

w 

are thus somewhat weak. It has been noted that the use of the leading 
term for ML^ and for (kML^ + ML^) is somewhat better than for kL^ 
which implies that the similarity parameter for d-type terms (see 
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eq 3.7) is of different form than for 8-type terms. 

3.3.1 Limiting Forms for A and r 

Series expansions for A and r in l/K give: 

A = - 1 - ~ I- 'VA. + O0-/K 2 ) 

1 + y 2(2 - l)/7 

r , M + 0U/K 2 ) 

2 + V2(7 - l)/7 

These are the expressions used in deriving (3*8) and are the limiting 
forms for K » 1 given by Chernyi (1961, p. 189) . The strong effect 
of 7 on A is apparent as 0.0557 < A < 1 for 5/3 < 7 < 1. 

3.3.2 Effect of Wave Reflections on Aerodynamic Forces 

For no wave reflections X = 0 in equations (3.5) such that: 


ML X = = 0 


kML 2 = kML^ = kMM^ = k 2 !^ = k 2 MM^ = FK t (3-9) 


These are the results given by linearized piston theory (see, for 
example, Ashley and Zartarian, 1956) multiplied by FK t which accounts 
for local external flow effects. Series expansion for FK for K » 1 

T 

gives. 


IK 


T 


,1 7..+.1L . K + - -l2_t.il ... 

V2(7 - l)/7 V 2 (7 - l)/7 


3 7 2 - 67 - 1 

7(7 - 1) (7 + l) 2 
— 


1 + 0 (l/K 5 ) 



Thus all aerodynamic terms are singular as 7 -» 1 and wave reflec- 
tions are essential for the theory to yield reasonable results as 
7 -» 1. This singular behavior has been discussed by Miles (i 960 ). 

For M0 -><*>, A -»0.l4 for 7 = 1.4. The modification of the 
w 

above results for kML^ resulting from the inclusion of wave reflec- 
tions is, from ( 5 * 5 ) > a multiplying factor of 1 - 2 \ + 0 (l 2 ) = 0.72 

for 7 = 1.4. Thus the correction to lift curve slope can be more than 

25 percent for large values of M9 for 7 = 1.4 and a "local flow 

w 

theory" which would neglect this effect could be considerably in error 

for M9 > 1. 
w 

Omitting the wave reflections also gives ML^ = MM-[ = 0 which is 

also the case with piston theory in which the local pressure depends 

only on the local effective piston velocity. These terms are essentially 

C and C , respectively, which result normally from downwash lags. 
i • SI* 

a a 

The wave reflection process is one mechanism for generating an effect 
like a downwash lag, resulting in nonzero values for these coefficients 
when the reflections are included. A further discussion of d effects 
at hypersonic speeds had been given by Ericsson ( 1968 ) . 

3.4 Comparison With Other Theories 
Chawla (1958) has given the corresponding aerodynamic forces for 
third order piston theory for a wedge as: 



k 

kMM^ = ^ kMLg 


(3.10) 
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For M9 « 1 expansion of equation (5-6) gives , 

w 


kMLg = 1 + ^ K 




+ 0(K 5 ) 


which agrees with (3-5) to 0(K). The coefficients of are slightly 

different, as they should be, as piston theory is based on the series 
for an expansion wave which differs slightly in the third order term 
from the series for a shock wave (Lighthill, 1953) • For small K the 
effects of wave reflections become small so that and also 

become negligibly small. Thus, it is demonstrated that piston theory 
is a small K subcase of equations (3*5) • 

Van Dyke (1954b) has shown that the Newtonian plus centrifugal 
force (sometimes called Newtonian snowplow theory) is one limit of the 
hypersonic small disturbance theory. For Newtonian conditions, 7 -» 1 
and K » 1, the leading term of equations (3-8) give the results for 
Newtonian theory with centrifugal forces included. Letting 7=1 and 
using equations (3*8) it can be shown that results agree with those of 
Aroesty, et al. (19 66) which were obtained by an expansion about 7=1 
in a different fashion. The results for the wedge (3*8) can be regarded 
as a generalization of the Newtonian plus centrifugal force theory to 
arbitrary 7 and somewhat small values of K. 

The pitch rate derivatives are often computed on tangent wedge 
basis, or Newtonian theory omitting centrifugal force effects, by 
assuming that the local pressure iB that of a wedge having the same 
slope as local effective slope of the airfoil including motion effects. 
The effective slope of a lower surface is given by: 
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9 = 0 + h + 0 + x0 

e w 

where superscript dot indicates — , and V is used as the reference 

dt 

external velocity consistent with small disturbance theory. Using 

',= (’* 1) 9 e 2 

and performing the appropriate integrations, the relations (3-9) are 
obtained with kMLg = (7+1) K. Thus the dependence upon 7 and K 
is considerably different than that given by (3*8) for the rate deriva- 
tives. The value of a tangent wedge method for large K for the 
dynamic derivatives would be heuristic in view of this comparison. 

Piston theory, however, is one form of a tangent wedge theory which is 
satisfactory for K « 1. 

3.5 Flutter of a Typical Section 
An analytical, model that has been extensively used to examine 
flutter trends is a rigid, two-dimensional airfoil restrained in an 
airstream by torsional and translation springs . The dynamical equations 
and an exact solution for the neutral stability boundary based on piston 
theory aerodynamics are well known (for example, Bisplinghoff and Ashley, 
1962, chapter 6) . No such exact solution can be obtained using aero- 
dynamic forces that vary with frequency such as equations (3-2) and the 
flutter condition must be sought by an iterative or other technique. 
However, if frequency independent aerodynamic forces such as (3*5) are 
used, the following exact solution results: 
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The above solution can be applied with aerodynamics from (5.5) > 
(5-8), or ( 5 . 9 ) . If piston theory (5*10) is used, then the same result 
as Chawla (1958) is obtained. The present solution ( 5 .11) includes 


ML^ and MM^, which are zero in piston theory, and thus retains terms 

of one higher order in k. The solution ( 5 .11) is a function of the 

% 2 

structural parameters — , x , x^, r^ and of the aerodynamic 

a 

similarity parameters K = MQ^, 7, and (iM only. The parameter pM 



2b 

__ o 

gives an equivalence between Mach number and altitude (p = m/p c ) which 
has previously appeared in piston theory results (see Bisplinghoff and 
Ashley 1962 , p. 253). 

Some results are presented in figure 3.2 to indicate the relative 

0 c 

importance of ML, or a terms, and kML^, or q = •=- terms. The 

effect of ML^ is apparent at low values of pM and as can be seen 

from (3-H) where it appears only in the ratio ML^/pM. However, 

kML^ has an effect even for large values of pM as it appears in (3*11) 

independent of pM. (Similar results have been presented for low speed 

flutter by Lamboume ( 1967 ) • ) The results are accentuated for large 

values of co, /<u and inclusion of kML.' is essential for flutter 
hj a 4 

prediction. A decrease in 7 to 1 (not shown) also affects the result 
by about 10 percent. It is also of some significance that at flutter, k 
decreases rapidly as pM Increases, and is small for large values of 
pM; the approximations ( 3 . 8 ) are thus applicable for practical cases of 
large pM. 

One consequence of the different effects of 7 on each of MkL^, 

MkLg, and ML^ as indicated by the similarity parameters of equations 

3.8 is that a combined similarity parameter for flutter is not apparent. 

Thus, boundaries such as given in figure 3*2 must be calculated for 

each value of 7 and M0 . 

w 

With an attached bow wave, the upper and lower surfaces of the wedge 

are independent. The preceding development can thus be applied for 

determining the effects of initial angle of attack for < 9^ and 

a « 1 by taking half of (3*2), (3*5)> or (3*8), using 9 =9 + a 

v c w u 
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for the lower surface, and adding similar results based on 0 c = 8^ - 
for the upper surface. For K » 1, the leading terms of (5*8) are 
proportional to K and the resulting expressions involve only the 
average K, which is and there is no effect of angle of attack on 

any of the aerodynamic forces. The effect of the terms O(^) would give 
a very small modification of this result as a Q « 1. From piston 
theory for M0^ « 1, including initial angle of attack gives (Chawla, 


1958): 


kMLg = 1 + 


(7 + 1 ) 


„ + K 2 . K>‘ 
K + 5 + —5— 


where |k ± Ma Q | < 1. For K -> 0 and Mo^ = 0.5, for example, the 
effect of Moq is 15 percent for y = 1.4. Thus the effects of a small 
initial angle of attack for a thin wedge airfoil varies from a modest 
effect for <1 to no effect for M0^ » 1. The initial flutter 
speed gradient can, however, be sizeable if M is very large and 0 

w 


is small. 
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4. LINEARIZED PERTURBATIONS ABOUT STEADY FLOW CONDITIONS 

A two-dimensional flat surface with a sharp leading edge and exposed 
to a supersonic freestream flow of a perfect gas at a mean inclination 
angle of 9 W is considered while undergoing a specified time-dependent 
motion (figure 4.1). Three types of simple-harmonic oscillations are 
considered separately — pitching about the leading edge 8, translation 
normal to the mean position of the surface hy , and translation parallel 
to the mean position of the surface h x . The motion is considered to be 
an infinitesimal perturbation about the mean or steady-flow condition 
which leads to a linearized analysis describing the region between the 
moving surface and shock-wave boundaries beginning at the wedge tip and 
extending a preassigned distance downstream. Linearization of the 
governing equations permits the application of the boundary conditions 
at the mean position of the moving boundaries and also permits treatment 
of any rigid-body, time-dependent motion of the flat surface by super- 
position of the three motions considered. There are essentially three 
criteria that govern the amplitude limits for "small" oscillations: 

(l) fraction of Aq = - 0 W fraction of 0 W (or for very small 

0 W , fraction of p,*,), or (3) fraction of 0^ - 8 W where 0^ is the 
detachment angle. It might be noted that only compression surfaces are 
treated; however, in a linearized analysis of an expansion surface, the 
linear theory of Garrick and Rubinow (1946) can be applied considering 
the mean surface flow to be the effective freestream. Also for 
application to wedges, the base pressure is assumed constant and thus 
does not enter into the motion-related force and moment coefficients. 



SHOCK 
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Figure 4.1.- Oscillating flat surface and coordinate system 
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The hypersonic small distrubance theory treated in Chapter 3 also 
considers perturbations about the steady flow for thin wedges and hyper- 
sonic Mach numbers. Here both supersonic and hypersonic Mach numbers 
are treated, both large and small inclination angles are considered, 
and the fore-and-aft translational motion h x , which is not significant 
for hypersonic small disturbance conditions, is also analyzed. 


4.1 Development of Perturbation Equations 
The governing equations written in divergence form for the x-y 
coordinate system (e.g., Leipmann and Roshko, 1957 > Chapter 7) : 


continuity: 


|| + + = o 


(4.1a) 


x-momentum: 


+ a( p ^_pu 2 ) , + _ o 


(4.1b) 


y-momentum: 


d(j5v) + d(p.v) + d(p +_pv 2 ) _ 0 


(4ilc) 


energy: 


a 

St 


£. + (/ ~ -U - p(u^ + v^) 
y 2y 


+ It [ u ( p + i2 2T il 5(u2 + 

+ Ijp + ^27^- + ' ?2 )) 


= 0 


(4. Id) 


where the barred quantities are dimensional. 
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Rmfl.1 1 variations about the steady flow are considered by substitut- 
ing in equation 4.1: 

a 

P = P 0 + « Pi U,y) e 

- - - / ico t 

p = p Q + £ p x (x,y) e 

- — , — / v iit 

u = Uq + e U X (x,y) e 

- — / s iai t 

v = e v 1 (x,y) e 

where € is a small quantity; is the frequency of oscillation; 
p-^, and so forth, are complex-valued perturbations in the flow variables 
(subscript 1); and p 0 , and so forth, pertain to the steady or mean 
wedge flow. Neglecting terms 0(e2) and higher gives, for the 
y-momentum equation, for example 

_ _ Sv Sp, _ _ 

p o “o S + S? + p 0 v i * 0 

The equations are simpler in form if normalized by the shock- layer 
variables p Q and Uq and if the shock- layer wavelength, u q /cd is used 
as the characteristic length. Setting p^_ = P-^/(pQ u^ 2 ), Pq =■ P±/ Pq> 
u-l = Uj/uq, v 1 = v-l/uq, x = xcB/uq, and y = yw/u Q 

the resulting equations are: 

y-momentum : 

dv dp 

S + ^ + 1 ^l’° (4 - 2a) 

x-momentum : 

du dp 

+ “5E + iu i = ° 


(4.2b) 
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energy: 


Sp to (M q Pl - u x ) 

^_ + _ _^. + 1 * = 0 


(4.2c) 


B„ 


continuity: 

dp^ du 


i ' & 'i 

sr + sr + y + 1 pj 


= 0 


(4. 2d) 


where B Q 2 = M Q 2 - 1. Equation (4.2c) is obtained by subtracting (4.2b) 
from the complete energy equation. Note that the continuity equation 
(4. 2d) is not needed except to calculate p and will be disregarded 
here. Thus, (4.2a-c) is a system of three simultaneous, linear partial 
differential equations for p^, u,, and as functions of x and y 

and with complex, constant coefficients. 

4.2 Boundary Conditions 

4.2.1 Boundary Conditions at the Wedge Tip 

For the assumed attached leading edge shock, the shock displacement 
boundary condition at the wedge tip for pitching motion is 

8 (0) = 0 (4.3a) 

for plunging motion normal to surface (omitting h^. e icu 4) 

8 (0) = - cos Aq (4.3"b) 

and for translation parallel to surface 

8 (0) = sin A Q (4.3c) 

4.2.2 Boundary Condition at the Wedge Surface 

The tangent flow condition at the wedge surface (subscript s) for 
infinitesimal pitching motion is (also omitting 9 e ico 1) 
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v i, s = - d + lx ) 

for plunging motion normal to the surface 

(4.4a) 

V l,s » - 1 

and for translation parallel to the surface 

( 4.4b ) 

V.. =0 

l,s 

(4.4c) 


4. 2. 3 Boundary Conditions at the Shock Wave 

The perturbations of the flow variables at the moving oblique shock 


wave are discussed 

in Section 8.2. At the shock boundary 

(subscript 6 ) 

the changes in the 

flow variables are given by equations 

( 8.2 - 8 . 5 ) in 

the following form 

*1,8 = P p » + 1 P v 6 

(4.5a) 


u l,8 

(4.5b) 


v l,8 * 13 + i V v S 

(4.5c) 

where 

3 . g oob lo 

(4.5d) 


and Pp = 8 and so forth, are derivatives of the shock relations 

evaluated at the steady or mean conditions. 


4.3 Approximate Solution by the Method of Integral Relations 
For the one-strip integral-relations approach used here, a linear 
y-variation of p, u, and v is assumed. For example 


p x (x,y) = P s (x) + [p 6 (x) - P s (x)] y/ 8 (x) total 


(b.6) 
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where subscripts s and 5 refer to the surface and to the shock, wave. 


respectively, §( x )-total = + ^( x )/ cos i s 4be total shock layer 

thickness, and p and p 5 are perturbation variables (subscript 1 
omitted). Equation (4.6) and similar relations for u and v are used 
in (4.2a-e) and integrated in the y-dlrection. Sample terms in (4.2) 
are: 



dy = p s (x) + P 5 (x)] S(x)/2 + 0(e 2 ) 



^P n (x,y) 

~V~ ^ 


P 5 (x) - p s (x) + 0(e2) 


and using Leibnitz's rule, 


f 5 d r & / s , d5 , , 

J — ay = / Pl (x, y ) ay - - p 8 (x) 


0 

dp (x) dp fi (x) 
s + o 


dx 


dx J 2 


5 0 (x) 


i r T d 5 n (x) o 

+ 5 [p s (x) - p 8 (x)] — 55T - + 0(e > 


This procedure leads to: 


d&. 


(p 6 + p s + u 6 + U s) - b (p 6 ' Ps + u 6 " u s } dT + i(u S + U s } = ° 


d5 n p 

s < v 5 + v s) - 6T (v 6 - v s) d3T + 6l (p 5 - p s) + ^8 + v s) = 0 
0 u 


d5, 


d (p* + Pj " (P fi - P B ) — + 
0 


0 . 2 


VP 6 T ±V " 8^ v *8 " *6' dx- ' ' 2 (v 5 " V s } 

5 0 B 0 


+ \ [ M 0 2 (p 6 + p s) - ^b + u s)] = 0 
Bo 



?4 

Now 

& 0 (x) = x tan ?\q 

d. 5 q (x) 

— aS ’ tan N) 

Furthermore, the 6-terms are given by the boundary conditions at the 
shock wave (eq. 4.5)* For example, 

P 8 (x) = P p p(x) + i P v 8(x) 

and thus 
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x 


B„ 


u + x P ^ 
s p dx 


+ 





+ i 


2 V 

V 

{ lx o 

u 

1 V — 

2 v s 

§ a ■ — 

2 2 v 

Bq tan Xq 

■4 " 1 X q 

4 B 0 ) 

— 1 x p 

B o _ 

B o 2 tan \ 


(4.7c) 


with 


cLS 

dx 


p/cos X Q 


(4.74) 


The partial differential equations (4.2) and their boundary conditions 
have been transformed to four inhomogeneous, complex, variable- 
coefficient, ordinary differential equations, in p s , u g , 6 and p. 

The known surface velocity, v s (x), appears as a forcing function, and 
the required initial conditions are also contained in (4.7) as will be 
subsequently shown. 

4.3.1 Initial Values at the Wedge Tip 

Taking the limit as x 0 in (4.7) gives 

v (0) V 

p(o) - -V- - i T- 5 (°) (4.8a) 

P (3 

p s (0) = P p p(0) + i P y 6(0) (4.8b) 


u s (0) = U p 3(0) + i U v 6(0) 


(4.8c) 
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These conditions are the exact quasi-static linearized conditions from 
oblique shock theory. 


4-. 3 . 2 Initial derivatives at the wedge tip 

Differentiating (4-. 7 ) with respect to x, taking the limit as 
x -* 0, and applying (4-. 8) gives an algebraic system which can be solved 
for 

d( 3 /dx|x = o> d P s /Hx=0> and du s /dx Uo 


The results are as follows : 


dx 


x=0 


- T<(1 - B 0 2 tan2 Xq) 


dv 

£ 

dx 


x=0 


- i tan Xq 
2 


B q 2 tan X Q 


b: 


(0) + i (m^ P v -U y + V y £ tan A Q ) B(0) 

( 4 . 9a> ) 


V B 2 P B 2 

+ I Mq*" P„ -u„+ rrrV + — V q tan A 0 + J~~\~ i P(°) 


3 sin 2 v p ''0 cos A, 


where 


A = v p + P p B o tan A 0 


(4.9b) 


and 


dx 


x=0 


-P 

p dx 


+ tan X, 


+ i 


x=0 
tan Xq 


dv 

E 

'0 dx 


x=0 


2 P 


, s (0) + i y v 6(0) + Vp+isi- P(0) 


(4.9c) 


du 
s 

dx 


dp 


x=0 


s 

dx 


x=0 


- i u s (0) 


()».9d) 
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4.3*3 Definition of Force and Moment Coefficients <• 

The motion-related aerodynamic force and moment coefficients are 
defined for a single (upper, see fig. 4.1) surface by 

— T c — - - — 2 

L 1 = - / p(x) dx = 4q ck 

J 0 


n 1 


2(L X + ±Z 2 ) -2- 


+ (L'^ + iL'^J 9+2 (Lj + iLg) 



(4.10a) 


re \ h 

' = / x p (x) dx = - 2q^ c 2 k 2 2 (M’-^ + iM’g) =?■ 

J 0 L 

h~ 

+ (M* 5 + iM'^) 6 + 2 (M' ? + iM'g) ~ 


lent 


(4.10b) 


These definitions are similar to those often used in two-dimensional 
supersonic flutter analysis (e.g., Garrick and Rubinow, 1946). Here, 
however, the force L is normal to the surface, rather than in the 
lift direction (indicated by ~ over the L) , and subscripts 7 an( i 8 
pertain to fore-and-aft translation parallel to the mean surface position. 
(The symbol N and the subscripts 5 and 6 are often used for unsteady 
flap hinge moments and for pertaining to a flap, respectively, see 
Garrick and Rubinow, 1946 . ) It might also be noted that no forces 
along the surface are generated, since the effects of viscous skin 
friction have been neglected and the base pressure has been assumed to 
remain constant. The aerodynamic coefficients for a wedge are given in 
terms of the above definitions in Section 8.3* 
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4. 3.4 Low-k Solution Based on Initial Conditions 

The surface pressure for small values of x (i.e., ^/uq) 
approximated as 


P s ( x) = p s (0) + 


x=0 


can be 

(4.11) 


Using the boundary conditions (4.3) in equations (4.8) and (4.9) and 
integrating for the forces and moment coefficients as defined by 
equation (4.6) gives for small values of k 




p 0 dp s ,hv 


dx 


x=0 


~ Pn U n 

k L 2 . * w s - - p (0) 


E , . „, 3 . . P J*L Ps;e ( 0 ) 


, 3 , P 0 U 0 ip s,0 

k - £ k M*4 — 


x=0 


r _ 3 M « ^2 dPs ' h * 

L 7 i M 7 "2 dx 


x=0 


k Lg = k M'g = 


P 0 U 0 


2 


(4.12) 


where the above coefficients are given for a single exposed surface. 
Evaluating equations (4.12) using equations (4.3, 4.4, 4.8 and 4.9) 
leads to the following algebraic expressions for the coefficients 
involving only the properties of the steady flow p Q , Uq, and so forth, 
and the shock- wave derivatives P^, V^, and so forth. 
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~ p o u o 

k L 2 = rvp (p p " A i cos 


k 2 Li = 


po u , 


0 0 


3 2 v p p 


k L 0 = 


P 0 U 0 


A, sin A- 


8 ~ 2 V ~1 0X11 A 0 


m Pn 4an A/-. r -i 

L 7 = | v ' p~A~" L (V v + Mo 2 P p sin \)> A 1 - P p A 2 sin *oj 


L-l = - cot Aq 1^ - p Q A^/2 


~ p o u o 

k H = 


[tan A 0 + P p (1 - B 0 2 tan 2 A^/A + aJ 


(4.13) 


where 


A = P Q V - P V. 
1 P v v P 


A_ = U. V - U V n 
2 p v v P 


a 3 - tan \>< 1 - vTa 


sin 


V - p P <v “o 2 v T e B ° 2 tan *°> 


and A is given by equation (4.9b). 


4 . 3.5 Numerical Solution of One-Strip Equations 

To determine the variation of the force and moment coefficients with 
k, equations (4.7) are numerically integrated using a computer sub- 
routine for first-order, real, ordinary differential equations. The 
complex equations (4.7) are converted to a system of eight coupled, 
first-order equations relating the real and imaginary parts of 



to 


5, p, u g , and p g . Initial values and derivatives at a small value of 
x are calculated from (4.1l), for example, which requires using (4.8) 
and (4.9) for starting the numerical integration of the differential 
equations. As k = ^ ~ = Uq x/ 2 , running integrals over x of the 

Vco 

real and imaginary parts of the pressure give the variation of the 
aerodynamic coefficients with k. Using the definitions (4.10) the 
coefficients are related to the integral of the pressure for a single 
exposed surface for the 

P r rx 


J 1 or 7 = ' J I 0 Be K or h* (l) ] 


as 


k L 


'2 or 8 


P 0 u a r x 
~2x 


~ J 0 


M. 


2 p 0 r x 


r- («)]«« 


1 7 x 3 J 0 i h - or 


k M, 


p_u„ r x 


2 or 8 


i dl 


(4.14) 


where | is the variable of integration in the x-direction. Similarly, 


* 2 % - - f2 sr / 0 Be [> <5>] « 


k H m -^r f 0 In [ p e (s) ] d{ 


for pitch, 



4l 


k2 ■s - - pj tJ 0 x K (t) ]‘ « 

kMJ = J Im j^p Q (6)J I dl (4.15) 

A computing program is given in Section 8.4 which calculates twice the 
above integrals. The trapezoidal rule is used for numerical integration 
of the surface pressures in the integrals to determine to coefficients. 

4.4 Results and Discussion 
4.4.1 Comments on the Method 

One result of applying the method of integral relations is a change 
in the region of influence of the governing differential equations. 
Disturbances propagate along the normal coordinate rather than along 
characteristics and the description of wave-type phenomena is thus 
altered. The application to steady supersonic flow over pointed bodies 
has been discussed by South and Newman, 1965* It was found that 
although the results of the method differed in detail, the exact wave 
behavior was approximated as a result of compensating effects. The 
integral method has also been applied to the perturbation of plane 
entropy layers which involves extensive wave effects with favorable 
results (George, 1987)* As discussed in McIntosh (1965a), and Hui 
(1969b), the acoustic waves generated by the surface and shock wave are 
important factors in determining the unsteady forces on the wedge. The 
description of this wave phenomena by the one- strip integral method 
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(eqs. 4 . 7 ) is investigated by applying it to a steady angular perturba- 
tion of a wedge surface. The results are compared in figure 4.2 with 
the exact solution from hypersonic small disturbance theory of McIntosh 
( 1965 a). The results of the one-strip method do not follow the exact 
wave pattern, but approximate it in a smoothed oscillatory fashion. 

This is, thus, an important limitation of the method if local details 
are important. However, it might also be noted that the example cited 
contains a slope-discontinuity. For cases with continuous slopes, this 
smoothing should be less severe. 

4.4.2 Low Reduced-Frequency Aerodynamic Forces 

4. 4. 2.1 Comparisons with other analytical results and experiment 

The exact solutions for low-frequency or stability-derivative-type 
of aerodynamics of Van Dyke (1953) and Hui ( 1969 a) are available for a 
symmetrical wedge pitching about an axis located on the chordline at any 
chordwise location. Direct algebraic comparisons with these results 
would be rather lengthy. The equations of Van Dyke (1953) were 
programmed and calculations made for an extensive range of M, 7 , and 
9 W . The results obtained were identical to the results from using 
equations (4.13) and the transformation equations of Section 8 . 3 . 
Furthermore, the results presented in the figures of Hui ( 1969 a) were 
also reproduced using the one-strip equations. However, no results for 
are available for comparison. Although not rigorously demonstrated, 
the results of the one-strip approximation (eqs. 4.13) are thought to be 
exact for the low-frequency aerodynamics. Further comparisons with the 
linear theory of Garrick and Rubinow (1946) and the second-order theory 
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Figure k.2.- Pressure perturbation downstream of a unit angle c han ge 
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of Van Dyke ( 1954a ) were also made by expanding the one-strip results 
in powers of 0 and obtaining identical results. Similar expansions 
for hypersonic small disturbance conditions agreed with the results of 
McIntosh ( 1965 a ,b ) . 

With the previously noted smoothing of the wave behavior by the 
method of integral relations, the degree of validity of these results 
may seem surprising. It was found by Hui ( 1969 a), however, that linear 
functions in x and y satisfied the perturbation equations for the 
pitching wedge and the caret wing. The linear y-function is the form of 
solution assumed by the one-strip method of integral relations (i.e., 
eq. 4.6) and the linear x- function assumed for the low frequency 
approximation (eq. 4.11). 

Calculated results using (4.13) are compared with the experimental 
data of Pugh and Woodgate ( 1963 ) for a symmetrical pitching wedge in 
figures 4.3 and 4.4. Generally good agreement is obtained with the best 
agreement for the thinner wedge and higher Mach number. 

4. 4. 2. 2 Selected Analytical Results 

The six normal-force coefficients (eqs. 4.12-4.13) for a single, 
inclined surface are presented in figures 4.5 and 4.6 for M = 2 and 
and for several values of the ratio of specific heats 7 . All coeffi- 
cients vary considerably with inclination angle, particularly for 
M = 00 , and approach infinity as detachment is approached. The coeffi- 
cients and k L'^ also change sign as detachment is approached. 

As the damping-in-pitch (about the leading edge) is related to k L'^ 
(i.e., see eq. 4.12), the damping-in-pitch thus changes to an unstable 
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value near detachment . This unstable pitch- damping has also been 
discussed by Hui ( 1969 a). 

The effect of y (figs. 4. 5-4. 6 ) is primarily to shift the detach- 
ment angle and thus has large effects near detachment even for small 
changes in the value of y. 

The fore- and aft-coefficients and k Lq are, of course, 

zero for 9 =0 and are somewhat small until the detachment angle is 

w 

approached. 

It appears that the detachment condition leads to singular behavior 

in each stability derivative. That such is the case for the static 

forces is apparent from the charts of Ames Research Staff (1955) which 

show that dp /d9„ is infinite at detachment. For sonic local embedded 

flow the static derivatives are not infinite, but are extremely large. 

For the small range of 0 for subsonic but attached flow, theories for 

w 

an infinite wedge are no longer valid as the subsonic flow field between 
the shock and the surface is not uniform for a finite wedge (e.g., 
Johnston, 1955 ). Some implications of the approach to detachment on 
the flutter characteristics of diamond airfoils have been discussed by 
Yates and Bennett (1971) • 

Generally, in linear aerodynamics, the presence of a singularity 
such as that shown here indicates a need for a more complete treatment 
of the nonlinear governing flow equations in order to describe a more 
realistic behavior. Here, however, at least up to conditions of sonic 
local flow, the exact infinitesimal amplitude results are obtained. 
Furthermore, the usual transonic refinement to linear theory retains a 
local x-derivative of the steady flow field as a variable coefficient in 
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the unsteady equations (e.g., Landahl, 1961 ); here the term is zero. 

The results presented here suggest that infinitesimal motion has limited 
practical significance near detachment since even small finite amplitudes 
result in rapidly varying aerodynamic forces. Near detachment the 
condition for infinitesimal motion is that the amplitude of motion must 
be small compared to 9^ - & v , where 0 ^ is the detachment angle. 

4.4.5 Variation of Aerodynamic Forces With Reduced Frequency 

4. 4. 3.1 Comparisons with other results 

There are no results known to the author for the variation of the 
aerodynamic forces with reduced frequency as computed from the super- 
sonic flow solution for the wedge of Carrier ( 1949b) or the extension of 
Van Dyke (1953). Furthermore, such computations would be a rather 
lengthy task. Consequently, typical results of the integration of 
equation (4.7) are compared for 9 W = 0 with supersonic linear theory 
for the flat plate (Garrick and Rubinow, 1946) are compared for small 

9 with second order theory of Van Dyke (1954a), and are compared with 
w 

supersonic small disturbance theory (McIntosh, 1965 a) in figures 
4. 7-4. 9 . In general, the one-strip integral method predicts the correct 
trands, but is accurate only when the frequency effects are small, such 
as in hypersonic small disturbance theory (fig. 4.9). Such a result 
might be anticipated in view of the previously discussed treatment of 
wave behavior by the integral method. 
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Figure 4.8.- Comparison of frequency effects for t/c = .1, M = 10/7 
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The results of integrating equations (4.7) for M = °°, y = 7/5> 
and 9 W = 25° are given in figure 4.10. The frequency effects are 
small for k « 1 which is the usual range of interest for hypersonic 
speeds . 

As previously discussed, the low frequency aerodynamic coefficients 
vary rapidly with inclination angle as detachment is approached. Corre- 
sponding frequency effects are presented in figure 4.11. A very strong 
frequency dependence even at low reduced frequencies is evident as 
detachment is approached. In view of the previous comparisons with 
other theories (figs. 4. 7-4. 9), these results may be considered only 
qualitative in nature. The trends indicated are important ones, however. 
For example, the coefficient k L'^ for 9 W = 22° (fig. 4.11) varies so 
rapidly as to suggest that the linear results may he physically unreal- 
istic. If a nonlinear analysis were required for this region, then the 
phenomena of subharmonics or higher harmonics may occur which are 
rejected by the linear theory. In fact the quasi-static results of 
Kuiken (1969) and Kacprzynski (1968) indicate that for some conditions, 
the unsteady pressure waveform is not a single harmonic of the frequency 
of oscillation, and shifts in the mean pressure level can occur even for 
moderate amplitudes of motion. Taking such nonlinear phenomena into 
account in a full dynamic analysis would require considerably more 
effort than would be required for a conventional linear analysis. A 
further discussion of these effects is given subsequently in Chapter 5. 







Figure k. 11.- Effects of k for several inclination angles, single 
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Rapidly varying coefficients may also lead to rapidly varying 
flutter speed-indices with parameters that affect the resulting k 
flutter such as the mass ratio (e.g., see Lambourne, 1967). 
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5. CALCULATION OF THE UNSTEADY FLOW FIELD BY A NUMERICAL 
FINITE DIFFERENCE TECHNIQUE 

5.1 General Considerations 

A numerical finite-difference calculation of the complete unsteady 
flow field of the wedge is presented in this chapter. The basic motiva- 
tion is to provide a means of calculating nonlinear aerodynamic forces 
to assess the range of validity of perturbation analyses and gain insight 
into the types of phenomena that may be encountered when the full 
nonlinear governing flow equations are treated. Of particular interest 
would be a method that could be applied to blunt bodies in order to 
evaluate the effects of oscillations in the subsonic embedded flow 
around the nose on the after body. The data of Muhlstein (1971) indicate 
that even low-amplitude, static perturbations may generate significant 
nonlinear effects in a transonic flow. Large frequency effects have 
also been demonstrated in Chapter 4 for wedge flows as the embedded 
flow field becomes transonic near detachment. 

As previously indicated in Section 2.2, several investigators 
(e.g., Barnwell, 1971 or Moretti and Abbett, 1966) have evaluated 
steady flow fields from time evolution of an assumed initial flow field. 
One might anticipate that suitable modifications for motion boundary 
conditions, unsteady flows could be treated. Moving boundaries, 
however, become quite difficult to treat numerically if the boundaries 
move with respect to the grid network. A key factor of the analysis 
used herein is a formulation of the problem in a body-fixed axis system. 
The grid network thus moves with the body and the body boundary condition 
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is considerably simplified. For the wedge in supersonic or hypersonic 
flow, only the near triangular region hounded by the vertex, by the 
moving shock wave, and by the wedge surface needs to be treated. Thus 
the flow field of interest is considerably more confined than a transonic 
or subsonic flow field. A body-axis formulation of the unsteady flow 
over a cone has previously been used in a perturbation analysis by 
Brong (1965). 

The use of a body- fixed axis system leads to the inclusion of body- 
force-type terms in the flow momentum equations that are proportional to 
0, h , 0, and so forth. Various types of motion can be treated by 

y 

performing calculations similar to the time asymtotic steady flow 
calculation. For example, pure pitch rate effects can be obtained by 
considering 0 = constant without applying the real world constraint of 


p t 

0 = / 0 dt. Similarly, if 0,0, and 0 are constrained appropriately 

J 0 

at each time step an oscillation can be executed. Oscillations, however, 

• •• 

would require a converged initial flow field for 0 and 0, and so 
forth, at t = 0. One of the chief merits of this procedure is the 
potential for determining steady flows aerodynamics, stability-derivative 
type aerodynamics, and oscillatory aerodynamics, all out of a single 
computer program, and thus potentially effecting considerable savings 
in the writing and debugging of programs. 


5.2 Governing Flow Equations 

The governing flow equations are written in divergence form in an 
(x 1 - y’ ) axis system that is fixed in the surface with its origin at 
the vertex. That is, the x’-y' axis system is the x-y axis system of 
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figure 4.1 except that it rotates and translates with the wedge. Here, 

primed quantities refer to the body- fixed axis system whereas unprimed 

quantities refer to the x-y system (fig. 4.1). Pitch 0 and both 

translational motions h x and h^ are treated, and are considered for 

combined motions. However, it is considered more natural for an aero- 

elastic problem to measure h x and hy relative to the undisturbed 

position as shown in figure 4.1, rather than in the body-axis system. 

The basic transformation equations for the momentum equations are 

discussed by Pai ( 1956 ) , and a somewhat similar development stated by 

Brong (1965). In this chapter quantities are nondimensionalized as 

x' =x'/c, y' = y' /c, t = t/(7 M /c) and u* = =- , v = =- , p = p/p ro , 

Voo »00 

and P = p/Cp^ V^ 2 ). 

The velocities in the axis system are related by 

u = U Q + u' + y* 0 

v = V Q + v' - x' 0 (5.1) 

where U Q , V Q are the velocities of the origin in the x-y system and 
are 

U Q = - (h x + cos 0 w ) cos 0 + (hy - sin 0 W ) sin 0 

V Q = - (h x + cos 0 W ) sin 0 - (hy - sin 0 W ) cos 0 (5«2) 


The substantial derivatives in the x* and y* directions are (Pai, 1956) 


du* du . 1 , 1 du . . du 

aF ’ Si * T0 + u ST + v W 


dv 1 _ dv 
dt ~ "St' 


u0 + u' 


dv 


+ u' 


dv 


(5.5) 
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Substituting (5*2) into (5-1) to obtain u and v and differentiating 
with respect to time and thence substituting into ( 5 * 5 ) gives the 
resulting momentum equations. The continuity equation in divergence 
form (Pai, 1956) is 


continuity 



(5.4a) 


Using (5.4a) the momentum equations can be written in divergence form as: 


x- momentum 


if ♦ s|r (P ♦ «»■*) ♦ 

- p cos 0 - hy sin 0 - y' 0 + x* 0 2 - 2 v' 9 J = 0 (5.4b) 


y- momentum: 


+ 'Sx 7 (p u ' v ') + (p + P v ' 2 ) 


-efe sin 0 + iiy cos 0 + x' 0 + y* 0 2 + 2 u' 0 ^ = 0 (5.4c) 


The transformation of the flow equations to arbitrary coordinate systems 
has been considered by Walkden (1966) who shows that the energy equation 
is essentially invariant under the transformation. For a perfect gas the 
energy equation written in divergence form is thus 


energy : 


It + p < u ' 2 + v ' 2) | + ^ 7 f u ’ p + ^ 27 ^ p ( u * 2 + v ' 2) 




2 7 

+ ^27 p ( u ' 2 + v ' 2 ) 


= 0 


(5-^d) 
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which readily reduces to the more familar form of (e.g., Sauerwine, 1965) 


£4 . 1 = o 

Dt p Dt 


( 5 - 4 e) 


5*3 Finite Difference Equations 

Motions of the wedge are treated by a numerical finite- difference 
calculation using the complete nonlinear flow equations in a body-fixed 
coordinate system. The explicit finite- difference scheme used is one of 
first order accuracy given by Lax ( 1954 ) for hyperbolic equations. It 
has been further discussed and applied to aerodynamic problems, for 
example, by deJarnett (1966), by Bohachevsky and Rubin (19 66), and by 
Bohachevsky and Mates (1966). The method is formulated by writing the 
governing flow equations in divergence form and replacing them with 
finite difference equations that contain implicit artificial viscosity 
terms. Although the region of physical interest is the nearly triangular 
shock layer beginning at the wedge tip and extending downstream, the 
hyperbolic stability criterion for the numerical solution requires that 
the flow field must be known at an initial line downstream of the tip. 
Furthermore, the technique does not consider a shock wave as a disconti- 
nuity but smears the discontinuity over a few mesh points both into the 
shock layer and into the freestream. This considerably simplifies the 
treatment of the moving shock boundary, but at the expense of flow field 
resolution. 

Subsequent to Lax ( 1954 ), many investigators have presented 
refined hyperbolic difference schemes, for example, Lax and Wendroff 
( 1964 ), Burstein ( 1964 , 1967 a,b), Moretti and Abbett (1966), Sauerwein 
(1966), Lapidus (1967), and Gourlay and Morris (1968). Burstein ( 1964 ) 
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and Emery (1968) have presented comparisons of results from several 
methods, and the hook by Richtmeyer and Morton (1967) and the proceedings 
of a recent symposium (NASA publ., 1970) discuss several methods at 
length. Variations of the Lax-Wendroff (1964) scheme such as that of 
MacCormack (1969) are often recommended for both temporal and spatial 
resolution over the method used herein. The use of Lax's method is 
considered warranted, however, as it is conceptually one of the simplest 
and is also used for an intermediate step in at least one of the higher- 
order methods (Bur stein, 1967b). 

The flow field is divided into a rectangular grid system for each 
time step as shown in figure 5«1 with constant values of At/Ax' and 
At/£y' as determined from the stability criterion as discussed later. 

The system of equations (5*1) to (5*4) can be written as: 

dA. 8B dC. 

-i + 3? + w + D i ‘ 0 (5 - 5) 

where 

A 1 = p 
A 2 = pu' 

A^ = pv* 


£ + pU 2 
y 27 



B 1 = pu« 
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B 2 = p + pu' 2 


B, = pu’v' 
3 


\ = p + W 2 ' y - 1 PU 2 u* 


C 1 = pv* 


C 2 = pu'v* 


= p + pv' 2 


= P + P u2 V ' 


D 1 = 0 


D 2 = - p h x cos 0 - iiy sin 0 - y' 0 + x' 0 2 - 2 v' 0 
D, = - p 2 u’ 0 + x' 0 + y* © 2 + ii sin 0 + ii cos 0 

3 [_ x y 


sin 0 + h cos 0 

y 


D 4 = 0 


,2 .2 . ,2 

r = u' + v T 
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The procedure of Lax (1954) is to replace time derivatives of time 


step k with the averaged forward difference , 


^l.m ~ [^ A lh+j,m + ^ A i^-l,m + ^H.m+l + ^Z.m-l] 


A t 


where k, l , and m are indices for t, x', and y' , respectively, and 
to replace (5 B^/Sx*)^ and (dC-^/dy 1 ^ with the symmetric difference 
expressions 

2 A x' 

and 

2Aj' 


and to replace 


1 

¥ 


(D i ) k with the averaged value 

(Dj )^ . + (Dj n + (Dj.)* .. 

v i'Z+l,m ' i'Z-l,m N 1 Z,m+1 


+ (D ± ) 


k 

Z,m-1 


Substitution of these expressions into (5-5) gives the finite difference 
algorithm. 


/ • \k+l 1 


(Aj - + (Aj)^ , + (A, ) k + (Aj)* 

' iy 2 +l,m ' 1 '2-1,111 1 Z,m +1 ' I'Zjin-l. 


A t 
2 A x' 


^ Bi ^Z+l,m " ^ Bl ^Z-l,m 


A t 
2 A y* 




A t 

“T 


L (D i ) m,m + ( D ih-l,m + (D i^,nH-l < D i>Z,m-l 

(i = 1, 2,3,4) 


(5.7) 
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Thus, at the k + 1 time step can he determined in terms of the 
known variables at the mesh points (k,Z+l,m), (k,2-l,m), (k,l,m+l), and 
(k,l,m-l) at the previous time step, k. This is a one step, explicit, 
alternating, difference scheme as the point (k,Z ,m) does not appear, and 
only values from the previous time step are used to compute the new 
point. With the determined, the fluid properties at a grid point 

can be readily determined as can be seen by inspecting the expressions 
for the A^, equations (5*6). 

Expansion of equation (5*7) "by Taylor's series about the point 
(k,l,m) gives: 



+ 


A x ' 2 
4 A t 



+ 


2 

A y ,4f 
FAt 



+ 


0 (A 2 ) 


(i 


1,2,3, 4) 
(5-8) 


The last two terms are somewhat analogous to terms that appear in the 
viscous momentum equations and are called "artificial viscosity" terms. 
These terms result in a smearing of rapidly varying quantities, such as 
fluid properties across a shock wave, over several mesh spaces, and in 
some instances excessive smearing or smoothing of transient phenomena 
can occur. Barnwell (1967) has applied a modification of this scheme 
by using the conventional Lax scheme (eq. 5.7) at odd time steps and 
using an alternate scheme containing no "artificial viscosity" terms 
for even time steps. The effective smearing of transient phenomena is 
thus reduced. In some of the higher order difference schemes. 
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artificial viscosity terms are often added for numerical stability 
(see, for example, Emery, 1968 ). 

5-3-1 Initial and Boundary Conditions 

The oscillating wedge in a supersonic flow with an attached shock 
is an initial and boundary value problem. The entire flow field over 
the grid network must be specified for t = 0. This is accomplished 
by initializing the flow field with one for a steady wedge at the equiva- 
lent angle of attack for the motion at t = 0. The difference technique 
is then applied until the flow variables converge holding the body 
forces constant (D^ in eq. ( 5 * 7 ) )• 

Each of the boundaries of the grid network requires special consid- 
eration at each time step, both during the initialization procedure and 
during computation of the unsteady flow field. As previously noted, the 
numerical technique does not consider a shock wave as a discontinuity 
but smears the shock over a few mesh spaces to satisfy the jump condi- 
tions. Sufficient freestream points must be retained such that the 
moving shock remains well within the grid network and such that the flow 
variables have approached freestream values. About eight additional 
grid points beyond the nominal location of the moving shock are required. 
The outermost points are held constant at freestream values of the flow 
variables . 

At the wedge surface the tangential flow condition in the body 
axis system is v* = 0. To calculate the body point at (k+1, l , l), 
an approximation is used that has been used in the past (i.e., Burstein, 
1964), referred to as the reflection approximation. An image point at 
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a distance Ay* inside the body is considered to have associated with 
it - (v'^ 2 ) k ' (p 1 , 2 ^ > ( u '2,2) k and A sli 8ht modification 

of this approximation is used here to allow for unsteady effects. At 
y' =0, the y- momentum equation (5.4c) reduces to 

•^7 = p h^ sin 0 + ty cos 0 + x* 0+2 u’ 0 J , y* = 0 (5-9) 

Thus a pressure gradient at the surface exists that is proportional to 
the motion (for the surface considered which has zero curvature). The 
pressure for the virtual point is thus modified by setting 

P?,-l - Pi, 2 - 2 Air' ‘(5.10) 


It is also noted that one consequence of the use of this approximation 
(e.g., deJarnette, 1 966) is 


dv 1 


0, y’ 


= 0 


The values for the initial line are calculated from the steady flow 
oblique shock theory based on the instantaneous angle of attack and Mach ■ 
number. This is a low reduced frequency approximation which requires 
that the reduced frequency based on the distance to the starting line 
be small, usually of order 0.1 unless frequency of amplitude effects 
are large. 

Insufficient data are available to calculate the last line of the 
grid network at each time step. Either an extremely large region at 
t = 0, or an extrapolation of the data is required. Here the last lin e 
is generated by an extrapolation of the at time step k to provide the 
data to calculate at k+1 from, 
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( 5 . 11 ) 


(A ± ) 


k 

2+l,m 


(Ai) 


k 

l ,m-l 


+ ( A i) 


k 

l ,m+l 


(Ai) 


k 

1-2, m 


Moretti (1968) has discussed the importance of accurately treating 
the boundaries in fluid flow problems and incorporates characteristic- 
type routines for treating the boundary conditions .for calculating 
steady flow fields in a time-asymptotic fashion (see, for example, 
Moretti and Abbett, 1966). The methods presented herein for treating 
the boundaries are approximate ones. This is considered an area for 
possible later development. 


5-3.2 Stability Criterion 

There is no general stability criterion for nonlinear hyperbolic 
equations in three independent variables. In numerical solutions such 
as the one being considered, a criterion based upon local linearization 
is used as a guide and if instabilities are observed the mesh is suitably 
adjusted to eliminate the instability. For linearized equations, Hahn 
(1958) has shown that the Courant-Friedrichs-Lewy (or CFL) condition to 
be a sufficient condition for stability for simplicial gird networks, 
that is, a grid network that would use only three base points in the 
previous time surface to determine the properties at mesh point for two- 
dimensional unsteady flow. The CFL condition states that the domain of 
dependence of the difference equations must include the domain of 
dependence of the differential equations, which is the forward Mach cone. 
Sauerwein and Sussman (1964) have discussed implementation of the CFL 
cirterion for various simplicial grid networks. 
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The grid network for equation (5-7) is nonsimplicial. A necessary- 
condition for the stability of linear, constant-coefficient equations is 
the vonNeumann condition which requires that the absolute values of the 
eigenvalues of the amplification matrix of the difference scheme must 
be less than or equal to one. Heie and Leigh (1965) have applied this 
criterion to several nonsimplicial networks and demonstrated that some 
nonsimplicial networks that satisfy the CFL condition do not satisfy 
the vonNeumann condition and would lead to instabilities in numerical 
calculations. Furthermore, Burnstein (1967a) has given an example of 
the calculation of a transient flow problem containing stagnation points 
and sonic lines, in which the vonNeumann condition was satisfied in a 
locally linear approximation, but instabilities were encountered in 
the calculation which used the complete nonlinear equations. In this 
case the nonlinear effects were considered to be essential to the 
instability, thus indicating possible limitations on a stability 
criterion based on linearization of the nonlinear equations being 
considered. Moretti (1968) has suggested that improper treatment of 
boundary conditions can lead to instability. Hence, the heuristic 
stability criterion used herein is the simpler CFL criterion with 
adjustments made in the grid network as the results dictate. 

The application of the CFL condition is illustrated in figure 5.2. 
To satisfy the CFL condition the diamond connecting the outer base 
points at time k must enclose base of the forward Mach cone about the 
particle path from point (k+l,Z,m). The circle of radius a At must be 
enclosed by the diamond formed by the outer base points. The ratio of 
step sizes is then for a square grid, 
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At 1 

Ax' - u* + v'| +~\J~2 a 


(5-12) 


where u' an d v* include the local translation of the body-fixed grid 
network relative to the freestream. The ratio of Ax' to Ay' is 
somewhat arbitrary. Here the mesh is chosen to be ’’square for simplicity 
and also for a slight reduction in the number of multiply operations. 

As a result of the stability criterion, selection of the position 
of the initial line and the grid spacing on it sizes the time step and 
grid spacing downstream. Consider the initial line to be located at 
Xq' and no grid points on the initial line within the shock layer. 

For a specified frequency of oscillation, the resulting value of k at 
Xq' is represented by k^ . The grid spacing is for the square grid, 

1C 1 
0 


Ax' = Ay' = 


x„* tan 7^ 


mjL - 1 


(5-13) 


and the corresponding value of k at the body point l is, 

k l = k xo + tan \) (5>l4) 

The total number of time steps for a full time step, N k> is determined 
from the period of oscillation, P, and the time step, 

P = 2 kt = N k At = 2rt (5-15) 


Substituting into (5-15) from (5.1*0 and (5-13) gives, 



(5-16) 
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For initialization or stability derivative-type aerodynamics, the 
number of time steps need only be sufficient for convergence of the flow 
field. The values of nq and xq' determine the value of (0c/2v), 
for example, at the body points downstream. 

5.4 Results and Discussion 

5.4.1 Discussion of the Numerical Finite Difference Results 

The effect of increasing time step size as fraction of the CFL 
criterion, on converged steady flow results is shown in figure 5*3 for 
0 W = 20°, M = 00 (m = 1000 used here) and y = 1.4. Stability was 
maintained up to slightly less the 0.8 of the CFL criterion, with the 
more accurate results obtained as the step size approaches the maximum. 
The results of figure 5.3 required approximately 250 time steps for 
convergence with somewhat fewer time steps required for the larger 
time step. Convergence proceeds essentially downstream from the initial 
line with points near the initial line coverging quite rapidly. The 
finite difference scheme used here is of first order and thus requires 
a fine grid network to maintain accuracy. In these calculations 16 
(unstaggered) grid points were used on the initial line inside the shock layer 
in order to maintain the limited accuracy of about ten percent. 

The results for several inclination angles are shown in figure 'y.k. 
The relative accuracy is about the same for each inclination angle. 

Thus the difference in pressure level with angle, which would be 
required for the static aerodynamic force derivatives, would have about 
the same relative accuracy. 
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The results for including hy motion are presented in figure 5*5* 
The results are comparable to the previous results. Inclusion of hy 

is essentially a change in wedge angle and Mach number and is thus 

• • 

essentially a steady flow calculation as hy does not appear in the 
momentum equations (5-4). 

The results of the finite difference calculation including 9 
motion is shown in figure 5*6 for M = 2 and », y = 1.4, and 9 W = 20°. 
For M = 2 an unstable damping in pitch was calculated by the perturba- 
tion method (figs. 4.5 and 5.11). The corresponding slopes of the 
pressure (fig. 5.6) is in the stable sense. However the pressure should 
extrapolate back to the static value for (9 c/2 V^) of zero which it 
does not in this case. This effect results possibly from a large value 
of x9 at the initial line in conjunction with the surface pressure 
being constrained to static values at the starting line. It appears 
that a preferable procedure would he to use the results of the perturba- 
tion method of chapter 4 for the starting line in this type of 
calculation. 

The results of the finite differaice calculation presented herein 
are shown to be of limited accuracy both from the use of a first order 
difference scheme and from relatively simplified treatments of the 
boundary conditions. The basic formulation of the problem appears to 
be a useable one, however, but the basic difference grid network is 
required to be so fine that a significant range of amplitudes, and so 
forth, cannot be treated in practical terms. For example, the number 
of steps required from equation ( 5 « 1 6 ) to execute a cycle of oscillation 
would be prohibitive for significant results. Second order difference 
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30 
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nite difference calculation 
motion 
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schemes such as that used by Barnwell (1971) can give good results with 
a relatively coarse grid network. 

For further development it is recommended that a second-order 
difference scheme be used. Explicit treatment of the shock boundary 
condition may also be necessary to reduce the numbter of grid prints and . 
increase the accuracy of the calculation. The region between the shock 
and body must thence be mapped such that the shock does not cut across 
the grid network and reduce the time step severely from the CFL criterion 
Furthermore, the treatment of the body boundary condition should be 
improved by using a characteristic -type method such as that used by 
Barnwell (1971) or Moretti and Abbett (1966). This would necessitate the 
derivation of the appropriate compatibility relations. The recent 
solution of Hui (1970) for the perturbation about the quasi-static flow 
condition could be used for the initial line to further improve the 
capability. 


5.4.2 Discussion of Quasi-Static Wedge Flows 

At high supersonic or hypersonic speeds, the value of k at 
flutter is generally very small, 0.01 or less. The resulting phase lags 
from motion are quite small. Although even the small phase lags may 
be important in the dynamics of an aeroelastic system, the general 
character of the flow may be assessed by considering instantaneous 
response of the flow field or a quasi-static flow. Consider a Taylor's 


series expansion of the steady wedge flow pressure coefficient 

dC B 2 C 2 d 5 C -3 C a k 

P p 8^ nfl> nfl 

= C + 0 + ir- + 

P 


c - = c p + 3e Ee+ — 2 - f- + — + — ?if + --(5.i7) 
p o ^W d0_ 2 5 b ^ 4 24 


w 


'W 


de 
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where the partial derivatives with respect to t, x, and y are assumed 
to he zero in accordance with the previous comments. Let 9 he described 
hy simple harmonic motion at infinitesmal frequency, 

9 = 0q sin tut (5*18) 


Substituting (5.18) into (5. IT) and using trigometric formulas for 
reducing sin 2 tut, and so forth, gives for the difference in pressure 
from the steady flow value, 


A CL = 


'8 2 G 


L^w 


8 4 C 0 n 2 

_£ + _^_£__ g + 


'8c 

3^ + 


W 

S 5 C 0 2 
P 0 
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cos 2 cut 



sin 3 “ft 


cos ^ cut + ... 


(5.19) 


Now from Ames Research Staff (1953) 

C P = 7TT [ sin2 p “ 1 /^» 2 ] (5 ’ 20) 

2 

where the relation between 3 and 0 W is a cubic equation in sin 3 
with coefficients as functions of 1^,7 and 9 W . The derivatives of 
(5.19) can be evaluated by repeated differentiation of (5. 20) and the 
cubic equation relating 3 and 0 . 
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Inspection of equation (5«19) indicates, first a zero shift if 
second or higher even-derivatives are large j and higher harmonic terms 
(i.e., other than involving 8 Cp/d 0 w ) occur if second or higher deriva- 
tives are significantly large. By evaluating these derivatives, the 
regions of M^, 7 , and 0 V where significant nonlinear effects of zero 
shifts and higher harmonics occur can he observed. 

In figure 5.7, C_ vs 0 is given for M = 2 and °° and for 

p W 

several values of 7. The derivatives through the fourth are presented 
in figure 5.8. (it might be noted that indeterminant froms of 0/0 are 
encountered as 0 W — » 0 and of increasing order for the higher deriva- 
tives. These were circumvented by letting 0 w e , a small value, for 
M = 2 and using C p = (7 + l) sin 2 0 W at 0 W = 0 for M = 00 ). The 
character of the derivatives (figure 5 - 8 ) is quite different for varying 
7 and M^. All derivatives become infinite as detachment is approached 
suggesting that nearing detachment leads to large nonlinear effects . 

The results of hypersonic small disturbance theory were given by 
Kuiken (1969) for M0 W = 2 and 5 which indicated large nonlinear effects. 
For small angles and high Mach numbers the nonlinear effects are also 
shown here to be large as the higher derivatives are larger in comparison 
to the lower derivatives. 

One implication of (5.19) is that the coefficient of the first 
harmonic contains contributions from the 3>5>**» derivatives. Compari- 
sons of filtered experimental data with a linear theory would be 
inappropriate if nonlinear effects are large. 
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It appears from these brief results that significant nonlinear 
effects are to be encountered near detachment at all Mach numbers and 
for large amplitude motions of wedges of small angle at hypersonic speeds. 
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6. CONCLUDING REMARKS 

To evaluate both the methods of treating unsteady flows for super- 
sonic flows and the resulting trends, several aspects of the oscillating 
wedge have been considered. The two perturbation methods-hypersonic 
small disturbance theory and the more general method-involve only a 
moderate level of effort for analysis and computation for the oscillating 
wedge. Further extensions to treat more general configurations would 
involve considerably more effort as the entire region of interest of the 
steady flow field must be known in detail. A nonuniform steady flow 
would lead to variable coefficients in the differential equations 
governing the perturbation quantities such that interpolation of the 
steady flow would be required for numerical integration. Methods are 
of course, available that describe the steady flow field in detail. It 
would appear that direct numerical solution of the unsteady flows may 
be practical for complex configurations along the lines of the formula- 
tion presented in chapter 5. However, further development is required 
to assess this possibility in practical terms. 

The trends presented for the wedge indicate that for inclination 
angles near detachment, large nonlinear effects may be anticipated at 
least for motions occurring at low reduced frequencies. The practical 
implications of detachment requires further investigation. In particular, 
the implications for three-dimensional wings and bodies from the two- 
dimensional wedge may be limited to large aspect ratio surfaces. 
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8 . APPENDIX 


8.1 List of Symbols 

A partial list of the principal symbols used with their principal, 
definitions is given here. Some of the symbols are defined locally 
where used. It should be noted that different nondimens ional factors 
are used in Chapter 5 than in Chapter 4. Some symbols also have differ- 
ent meanings in different contexts. 

A^,B^,C^,D^ see equation 5.6 

a local speed of sound 



b airfoil semi chord 

C pressure coefficient, (p - P )/ ( 3 00 

ir U 00 

c reference length 

F see equation 3-4 

h x displacement parallel to surface (fig. 4.1) 

hy plunging displacement normal to surface (fig. 4.1) 

i V-l" 


K 

K 

T 

k 

k,l,m 


L. 

J 



M0 

w 

Mt, where T(also p) is shock wave angle 
reduced frequency, c 

finite difference grid indices for t,x,y respectively 
(fig. 5.1) 


coefficient in normal force expression (see eq. 4.10a) 
i = 1,2,3, 4,7,8 

Mach number, v/a 

coefficient in moment expression (see eq. 4.10b) 

j = 1,2, 3, 4, 7, 8 
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m 


m. 


P v' P |3 

P 

q. 

r ,r q 

v 7 3 

r 

a 

t 

t/c 

U > U A 
V p 

u 

V 


V V p 


v 

X 


n 


w 


€ 

r 


airfoil mass per unit span 

number of grid points on initial line within shock layer 

^P s /^ v n> ^Pg/dp 
pressure 

dynamic pressure 
dp 5 /SV n , 3 p s /S3 

section radius of gyration about x^, units of b 
time 

thickness /chord ratio 
dU 5 /dV n , dU 6 /c)0 

velocity component in x-direction 
total speed 

aV 6^V aV 6^ P 

velocity component in y-direction 
coordinate distance (fig. 4.1) 

pitching axis, fraction of c" measured from leading edge 
coordinate distance (fig. 4.1) 
shock angle (fig. 4.1) 

2k(r n - l) where k is reduced frequency, n = 1,2... 
quantity defined by equation 4.9b 

shock displacement perturbation measured normal to steady 
shock 

shock layer thickness in y-direction (fig. 4.1) 

wedge angle (5 = l/2 tan -1 (t/c) for symmetrical wedge 

airfoil) 

perturbation parameter 
see equation 3*4 
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A 

*0 

co 

“h 

CO 

a 

P 

9 

0 

w 

M. 


ratio of specific heats 
see equation 3.4 

included angle between shock and surface (fig. 4.1) 

circular frequency of oscillation 

natural frequency in plunge 

natural frequency in pitch 

density 

pitch angle 

surface inclination angle 
mass ratio, m/p c^ 


Subscripts : 


f pertaining to flutter 

n normal to shock wave 

s evaluated at surface 

u,l pertaining to upper or lower surface respectively, for 

symmetrical wedge 

v pertaining to normal velocity of shock wave 

w pertaining to symmetrical wedge 

0 pertaining to change in shock wave slope 

6 evaluated at shock wave 

0 steady flow value in shock layer 

1 perturbation quantity 

0 ° freestream value 

Superscripts : 

' pertaining to pitch (or pitching moment) about leading 

edge or pertaining to body-fixed axis system 

~ pertaining to forces normal to surface 



Bar over variable indicates dimensional quantity 
Dot over variable indicates d/dt 
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8 .2 Perturbations of Flow Variables at a Moving Oblique Shock Wave 

The incremental values of p, p, u, and v resulting from perturba- 
tions of slope and normal velocity of the moving oblique shock wave are 
required as the boundary conditions for a perturbation analysis. The 
shock wave perturbation is measured normal to the steady shock wave 
position as shown in figure 4.1. The appropriate perturbations of the 
Rankine-Hugoniot conditions in a coordinate system aligned with the 
steady shock wave have been developed by Carrier (1949a) • They have also 
been given in the surface coordinate system used here (fig. 4.1) 
specialized to 7 = 7/5 (Carrier, 1949a and Van I)yke, 1953)* and for 
general values of 7 but in somewhat different nondimens ional units by 
Carrier (1949a) and Hui ( 1969 a) . 

With the normal velocity of the shock wave given by: 

V n = 15 ^ ( 8 . 1 ) 

J yjyl. 

the perturbations are written in the following form (omitting e ) : 


where 


p 5 = P p (3 + iP y 6 
p 5 = V + 1R v 5 

( 8 . 2 ) 

u B = y + iu v 5 

v 5 " V + 1V v 8 

P p = P y = 8p 6 /8v n , etc. 



evaluated at the steady shock wave position and subscript 6 denotes 
conditions just aft of the shock wave. Letting 


M n = M oo sin P 0 

2(M^ + 1) 

U = — p 

(7 + l)lt 

then the shock derivatives can be written as 


P = 


R = 


U = 


(7 + 1 ) PqUq Sln P 0 


4 p o 

^ + ^ M* sin p. 


“ U n Sin \) 


(8.4) 


V = U cos 
v n X) 


and 


P. = 


P 

— cos p_ 
u 0 0 


Ro = 


U„ = 


R cos p_ 
v 0 

U 

— sin 9 - p n P cos \ 

u Q w 0 v 0 


U 

V„ = — cos 9 - p. P sin 7- 

P u Q w 0 v 0 


(8.5) 


The above relations .reduce to those given by McIntosh (1965a) in the 
hypersonic small disturbance limit. It can also be shown with a similar 
development, that perturbations of the shock wave tangential to the 
steady shock position results in second order perturbations in the flow 


variables . 
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The steady flow values of pressure, density, and velocities within 
the shock layer are given in Ames Research Staff (1953) as functions of 

7, M m , and (3^ . The corresponding relation between p^ and is 

2 

also given as a cubic in sin p^ with coefficients as functions of 7 
and 9 . These relations are solved numeric ally for the steady flow 
field parameters. 

8.3 Calculation of Aerodynamic Coefficients for a Wedge 
from Surface Coefficients 

The unsteady aerodynamic coefficients for a single isolated surface 
have been derived in the body of this report (Chapter 4) for rigid body 
pitch and translation perpendicular and parallel to the surface. The 
coefficients and k were based on c, the length of the surface and the 
pitch axis was located at the leading edge. The coefficients for a 
symmetrical wedge at an angle of attack are given here in terms of the 
surface coefficients. The pitch axis is assumed to be on the wedge 
midplane or chordline. The derivation of these relations are briefly 
outlined and the results summarized. 

8 . 3 .I Force and Motion Transfer 

The forces and moments for a symmetrical wedge (subscript w) are 
related to the previously given coefficients for a single surface by 

/v 

L = ± L cos 6 
w w 

C = - L sin 6 
w w 

M' = ± M' 
w 


( 8 . 6 ) 



io6 


where the upper sign is for an upper surface. The pitch axis is dis- 
placed from the surface and the chords related by 


c = c /cos B 
w' w 


x A = x- cos 5 
0 0 ,w w 


y~ = - x. sin B 
0 0 ,w w 


(8.7) 


The perturbation displacements are related by 


h = ± h cos B - h sin B - x^ cos 8 9 

y y w w x w w °> w w 


H = h cos 8 ± h sin 8 - x^ sin B 9 

x x__ w y__ w 0 ,w w 


( 8 . 8 ) 


w 


w 


9 = ±9 

u, l 


where again the upper sign pertains to the upper surface. 

The force and moment coefficients are defined by equation 4- .10 for 
a single surface. By defining the coefficients in manner similar to 
(4-. 10 ) for the wedge, the coefficients for the surfaces and the wedge 
can be related through the use of ( 8 . 6 ) for the dimensional forces. 

Then substituting (8.7) and ( 8 . 8 ) into the expression for the surface 
coefficients, the wedge coefficients can be equated to the equivalent 
combination of surface coefficients . 

8 . 3.2 Coefficients for Wedge 

The results for the coefficients of lift for the wedge are: 


L i " (E iu + + (E 7u + E 7z } tan 5 w 

L 2 = (E 2u + E 2 1 ] + (E 8u + E 8Z ) tan S w 
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L 3“ (£ 3u + E 3z )/cos2 8 v 

K=< e 4u + l>W)/<=°= 2 

L 7 = (J 7u - SfP - - \P tan 5, 

^ - (*8u - *8 1> - (E 2U - K 2 1> *“ 5 * (8 - 9) 

and 

Mj_ = t(M[ u + + (Mj u + H^) tan 5 w ]/cos 2 & w 

*Z " ^ + *2l> + (M 8u + *S0 tan S w ]/cos2 6 w 

= (M^ u + M^)/cos 4 8 w 

M4 = (M^ + M^p/cos 4 S w 

Mj = [(Mj u - M^j) - (Mj_ u - M£j) tan 6j/ cos 2 & w 

%= C(M^ U - % z ) - (M^ - M^) tan 5 w ]/cos 2 & w (8.10) 


where the subscript has been omitted from the left hand side of ( 8 . 9 ) 

and (8.10). The relationship for reduced frequency for the wedge is 

k = k 7 cos 5 (8.11) 

u or l w v ' 

It may also be noted that for zero angle of attack, the coefficients for 

the upper and lower surfaces are equal and thus the coefficients 

Ly, kLg, My , and kMg are zero. Such would not be the case for nonzero 

angles of attack, however. 

For pitch axis locations other than x n = 0, the transfer rela- 

u,v 

tions are found from the same procedure to be the conventional ones 
(e.g. Garrick and Rubinow, 19^+6) and are 
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L l= 4 


and 


L 2 = L 2 


L 3 = L 3 - "Vl 
L 4 = L i _ 2 x o L 2 


L T = L 7 


l 8 " l 8 


\ " 2x o L i 

Mg = - 2 x q L 2 

M^ = - 2x q (M^ + Lj 

\ = % ~ % (Mg + L( 

^ = “7 " 


2x o L i) 


( 8 . 12 ) 


where X- = x rt in (8.12). 
0 0,w ' 


8.4 Description of Computer Programs 
The two principal programs used to generate the results presented 
in chapters 4 and 5. The FORTRAN computer program for numerically- 
integrating the complex system of equations (4.7) to obtain the aero- 
dynamic coefficients as a function of k from the perturbation is 
presented first. A library routine, INT2A, for integrating a system of 
real, first-order, differential equations is used to integrate (4.7) 
after expanding into real form. The pitch and two translational motions 
are treated separately. The program for finite difference calculation 
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of the unsteady flow field is then presented which can treat steady, 
oscillatory, or quasi-static motions. Pitch and the two translational 
motions may be treated separately or combined. 

These programs were written by the author for use on the Control 
Data 6000 - series machines at Langley Research Center using the RUN 
compiler and the Langley Research Center version of the SCOPE 3*0 
operating system. Approximately l4 significant figures were used in the 
computations. The compiler used permits the use of multiple arithmetic 
statement on one card when separated by the character $ . It might 
also be noted that the quantity 177700000000000000000g signifies an 
indefinite (undefined) quantity. 

8.4.1 Program for Perturbation Analysis 

8. 4. 1.1 Input 

Each case consists of a single card (80 characters) of identifica- 
tion for labeling the printout only, and list of variables in a NAMELIST 
called INPDATA. The FORTRAN variables and their definitions are as 
follows . 

FORTRAN VARIABLE Definition 

XM 

00 

THWD 9 , degrees 

w 

G 7 

XO starting value of k for beginning 

numerical integration 

XSTOP stopping value of k 

Cl k - increment for numeric all y inte- 

grating the differential equation 

k - increment for printing results 


SPEC 



3 
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8. 4. 1.2 Output 

The program lists each aerodynamic coefficient for each value of 
k requested (by SPEC) and the real and imaginary parts of the perturba- 
tion variables. The coefficients are also written on tape in coded (BCD) 
form for subsequent use by plotting program. 



o no 


in 


8.4. 1.3 Listing of Perturbation Program . 


OVERLAY!PF1STRP,0,0 ) 

PROGRAM PF1STRP( I NPUT=1, OUTPUT, TAPE7,TAPE5= INPUT) 

*$***** ijc********************** ***$******** **************** *********** 


* PROGRAM PF1STRP CALCULATES THE UNSTEADY FORCES ON AN OSCILLATING* * 

* INCLINED FLAT SURFACE IN A SUPERSONIC PERFECT GAS FLOW. * 

* THE LINEARIZED PERTURBATIONS ABOUT THE MEAN STEADY FLOW ARE * 

* CONSIDERED, WITH SOLUTION GIVEN BY A ONE-STRIP INTEGRAL METHOD. * 

* TWICE THE FORCES AND MOMENTS FOR A SINGLE SURFACE ARE PRINTED. * 


************************************************************************* 
COMMON /BLOC K1/CUVAR<9) ,0ER(9) ,VAR(9) , PDB , RODB , UDB , VDB , XM2 S , X BOS , 

4- RBOS»CVSR»CPLGI»CPTHI,CTNL2»DVSDXR»DVSDXI»RVDB,POV*UDV,VDV,RCNL 
+,DRO,DTO 

DIMENSION IDFNTI8) ,OATE (21 ,ELE1 ( 8) ,ELS2(8) , ERRVAL (8 I 
EXTERNAL DERSUB.CHSUB 

NAMELlST/INPCATA/XM,THWD,G,XCK,XSTOPK ,C IK, SPECK 
C 

110 FORMAT! 5E 1 6 • 8 ) 

109 FORMAT ( IH1 //* RIGID BODY FORE AND AFT TR ANSL AT I ON* ) 

108 FORMAT! /7X*K*AX*KSQL3P*6X*KLAP*AX*KSQM3P*6X*KMAP*9X*BR*8X*BI* 

4- 8X*PR*3X*P I*8X*UR*8X*UI*8X*DR*SX*DI*/) 

107 FORMAT! 1H1//* RIGID BODY PITCH*) 

106 FORMAT! /7 X* K*7X*L1 *7 X*KL 2*7X *M1 P*6X*KM2 P*8 X*BR*8 X*B I * 

+ 8X*PR*BX*PI*8X*UR*8X*UI*8X*DR*8X*DI*/> 

105 FORMAT! 1H1//* RIGID BODY PLUNGE NORMAL TO SURFACE*) 

10 A FORMAT! X F6 . A , AG 1 1 . A , 8G1 0 . 3 ) 

103 FORMAT! /7X*K*7X*L7*7X*KL8*7X*M7P*6X*KM3P*8X*BR*8X*BI* 

4- 8X*PR*8X*PI*8X*UR*8X*UI *8 X*DR*8X*D I * / ) 

102 FORMAT!//* P0B=*G16.8,* R OD B= *G 16 . 8 , * UDB=*G16.8,* VDB=*G16.8 
4-/* PDV = *G16.8,* R0DV = *G16 .8,* UDV=*G16.8,* VDV=*G16.8/) 

101 FORMAT! /2X*M2=*G16. 8,* P21=*G15.8,* R21=*G16.8,* V21=*G16.8,* B0=* 
4-G16.8) 

DATA FMTl,FMT2,FMT3/6H( 8A10) ,10H( 1H110A10) , 10H( * IOB*IA)/ 

VAR! 1, ..,9)=X,BR,BI,PR,PI,UR,UI ,DR,DI. RF=.5*V21*X 


RFWIND 7 

1 READ FMTl.IDENT $ I F ! EOF , 5 ) 999, 2 

2 CALL DAYTIM(DATE) $ PRINT FMT2 , I DENT , DAT E $ READ INPDATA 
PRINT INPDATA $ THW=THWD/57. 2957795 130823 

CALL WEDGE! XM,G»THW,XM2,P21,R2l»V21,B0, IEOBS) 

IF! IEOBS .EQ.O )G0 TO 3 * PRINT FMT3, IEOBS $ GO TO 1 

3 PRINT 10l»XM2,P21*R2l»V21,B0 

CALL SHKDERV! XM , G, THW , BO ,R21 , V2 1 , PDV , RODV, UDV , VDV , P DB , RODB , UDB , 
4- VDB) $ XM2 S= XM2**2 $ XB0S=XM2S-1. 

PRINT 102, PDB, R0C8, UDB, VDB, PDV, RODV, UDV, VDV $ PRINT 105 
R BOS= 1 . /XBOS $ RVDR=1 • /VDB $ PCON=-R 2 1 *V2 1 S C0NM=PC0N/3. 
C0NTM=-2.*R21/3. $ X0=2 . *X0K /V2 1 $ XST0P=2 . *X S TOPK/ V2 1 
CI=2.*CIK/V21 $ SPEC=2.*SPECK/V21 
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RIGID BGDY PLUNGE 

CVSR=0. i CPLGI=-1. $ CPTHr=0. $ 11=0 $ DR0=-CQS ( BO-THW) $ DI0=0. 
CALL INIC!DPDXR»DPDXI*XO, THW »B0 » PROO, P I 00 » T NL ) $ CTNL2=2./TNL 
CL10=-R21*DPDXR $ CMlP0=4.*CL10/3. $ CKL20=PC0N*PI00 $ XK=0. 
WRITE(7,U0)XK,CL10,CKL20,CM1P0,CKL20 
PRINT 106 SPRINT 104, XK, CLIO, CKL 20 *CM1P0»CKL20 

CALL INT2AI 1 I ,8,DUM,CI , SPEC , DUM , DUM , VAR , CUV AR , DER , ELE 1 ,ELE2 , DUM, 

+ ERRVAL, DERSUB.CHSUB, DUM) S XOO=VAR(1) 

C Ll=-R21* < PR OO+VAR ( 4 ) l/XOO $ CKL2= . 5*PC0N*{ PI 00+V AR < 5 ) ) 
CMIP=CONTM*(PRCO+2.*VAR(4) l/XOO $ CKM2 P=CONM* ( P I 00+2 . *VAR I 5 )) 

X SA V= VAR ( 1 ) $ PR SAV= VAR ( 4 I $ PISAV=VAR!5) $ XK= . 5 *V2 l*VAR (1) 

WRITE (7, l 10) XK,CL1,CKL2,CMIP,CKM2P 

PRINT 104, XK, CL1,CKL2,CM1P,CKM2P, (VAR( 11,1=2,9) $ 11=1 

7 CALL INT2AI 1 1 , 8 , DUM , C I , S PEC , DUM , DUM , VAR , CUV AR , DER , ELE 1 , EL62 , DUM, 

+ ERRVAL ,DERSUB,CHSUB,DUM) 

X S= VAR ( 1 ) S R XSO= 1 . / X S* *2 S DX = XS-XSAV 
R X = XSAV/ XS S RX2=RX*RX $ RX3=RX2*RX 
CL1 =RX2*CL1-R21*DX*<PRSAV+VAR(4) )*RXSQ 
CKL2=RX*CKL2 + .5*PC0N*DX ; <‘!PISAV + VAR( 5) )/XS 
TX1=VAR< l)+2.*XSAV S TX2=2 .* VAR ( 1 )+ XSAV 
CKM2P=RX2*CKM2P+C0NM*RXSQ*DX*(TX1*P I SAV *TX2 *V AR ( 5 ) ) 

CMlP = RX3*CMlP+CCNTM*RXSQ*DX*!TXl*PRSAV-»-TX2*VAR!4) )/XS 
X SA V=VAR ( 1 ) * PR SAV = VAR (4 ) $ PISAV = VAR(5) S XK= . 5*V2 1 *VAR ( 1 ) 

WRITE (7,1 10 )XK,CL1,CKL2,CM1P,CKM2P 

PRINT 104,XK,CL1,CKL2,CMIP,CKM2P,(VAR( 11,1-2,9) 

IF( XS.LT.XST0P1G0 TO 7 $ ENDFILE 7 $ PRINT 109 

RIGID BOOY FORE AND AFT TRANSLATION 

CVSP=CPLGI=CPTHI=DIO=0. $ 11=0 $ DRD=S IN ( BO-THW ) 

CALL INICIOPOXR.DPOXI , XO, THW ,BO,PROO,PIOO,TNL) $ PRINT 103 
CL70=-R21*DPDXR $ CM7P0=4.*CL70/3. $ CKL80= PC0N*PI00 S XK=0. 

WRITE (7, 110 )XK,CL70,CKL80,CM7P0,CKL8D 
PRINT 104,XK,CL70,CKL80,CM7PO,CKL80 

CALL INT2AI I I , 3 , DUM , C I , SPEC , DUM , DUM , V AR , CUV AP , DER , ELE 1 , ELE2 , DUM, 

+ ERRVAL, OERSUB,CHSUB, DUM) $ X00=VAR(1) 

CL7=-R2 1 *( P ROO+V AR ( 4 ) J/XOO $ CKL8= . 5*PC0N* ( PI 00+ VAR ( 5 ) ) 
CM7P=CQNTM*<PRC0+2.*VAP.<4) )/X00$ CKM8P=C0NM*( P 100+2 . *VAR ( 5 ) ) 

X SA V= VAR ( 1 ) $ PR SAV = VAR ( 4 ) $ PISAV=VAR(5) S XK=. . 5*V2 1 *VAR ( l ) 

WRITE (7, 110) XK,CL7,CKL8,CM7P,CKM8P 

PRINT 104, XK,CL7,CKL8,CM7P,CKM8P, (VAR( 11,1=2,9) S 11=1 

8 CALL INT2 A ( 1 1 , 8 , DUM , C I , S PE C , DUM , DUM , VAR , CUV AR , DER , ELE l , ELE2 , DUM, 

+ ERRVAL, DERSUB.CHSUB, DUM) 

X S = V AR ( 1 ) $ RXSQ=l./XS**2 * DX=XS-XSAV 

RX=XSAV/XS $ R X2 = RX*R X $ RX3=RX2*RX 

CL7=R X2*CL7-R21*DX* (PRSAV+VARI4) >*RXSQ 

CKL8=RX*CKL8+.5*PCON*DX*(PISAV+VAR( 5) )/XS 

TX1=VAR( 1)+2.*XSAV $ TX 2 = 2 . * VAR ( 1 ) + XS AV 

CKM8P=RX2*CKM8P+CCNM*RXSQ*DX*(TX1*PISAV+TX2*VAR (5) ) 

CM7P=RX3*CM7P+C0NTM*RXSQ*DX*(TX1*PRS4V+TX2*VAR (4) ) /XS 

X SA V= VA R ( 1 ) $ PR SAV = VAR 1 4 ) $ PISAV = VAR(5) $ X K= . 5 *V2 1 *V AR < 1 > 

WRITE!?, 110)XK,CL7,CKL8 , CM7P , CKM8P 

PRINT 104,XK,CL7,CKL8,CM7P,CKM8P,<VAR( I) ,1=2,9) 

IF! XS.LT. XSTOP)GO TO 8 $ ENDFILE 7 $ PRINT 107 
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RIGID BODY PITCH 

CVSR=-1. $ CPLGI=0. $ C PTH I=-l . $ 11 = 0 $ DR0=DI0=0. 

CALL INIC(DP0XR»0P0XI,X0, THW ,B0*PROO*PI 00* T NL I $ PRINT 108 
CKSL3P0=V21*PC0N*PR00 $CKL 4P0=PC0N*DPDX I $C KM4P0=4 . *C KL4P0/3 . 
XK=C. * WRITE(7,110)XK,CKSL3FO,CKL4PO,CKSL3PO,CKM4PO 
PRINT 10 A, XK,CKSI.3P0,CKL4P0,CKSL3P0,CKM4P0 

CALL INT 2 A ( 1 1 , 8 , DUM , C I , SPE C , DUM , DUM , V AR , CUV AR , DER , ELE 1 , ELE2 , DUM 
* ERRVAL,DERSUB,CHSUB,DUM) $ XOO=VAR(1) 
CKSL3P=.5*PC0N*V21*(PR0C+VARC4>) $ CKL4P=PC ON* ( P IOO+V AR ( 5 I) / XOO 
CKSM3P=C0NM*V21*< PR00*2.*VAR (4) > 
CKM4P=2.*C0NM*(PI00+2.*VARI5))/X00 

X SA V=VAR ( 1 ) $ PRSAV = VAR (4 ) t PISAV=VAR(5) * XK= . 5*V21 *VAR ( 1 ) 
WRITE(7*110)XK*CKSL3P»CKL4P,CKSM3P* CKM4P 

PRINT 104, XK, CKSL3P. CKL4P, CKSM3P, CKM4P, (VAR(I) , 1=2,9) t 11=1 
9 CALL INT2A ( 1 1 , 8 , DUM , C T , SPEC , DUM , DUM , V AR , CUV AR , DE R , EL E 1 , EL E2 , DUM 
+ ERRVAL ,DERSUB,CHSUB,DUM) 

X S= VAR ( II % R XSQ=1 . /X S**2 $ CX=XS-XSAV 
RX=XSAV/XS t RX2=RX*RX $ RX3=RX2*RX 
CKSL3P=RX*CKSL3P + .5*DX*V21 *PCON* ( PRSAV-*-VAR(4) ) / XS 
CKL4P=RX2*CKL4P +D X*R X SQ*PCON* ( P I S AV+ VAR C 5 I ) 

TX1=VAR(1I4-2.*XSAV t TX2=2.*VAR< 1H-XSAV 

CKSM3P=RX2*CKSM3P-*-C0NM*V2i*RXSQ*DX*{TXl*PRSAV+TX2*VAR (M ) 

CKM4P = RX3*CKM4P *2 . *CONM*RXSQ*DX* < TX 1*P I S AV+T X2*V AR ( 5 » )/XS 
X SA V= VAR ( I ) $ PRSAV = VAR (41 t PISAV = VAR(5) $ XK= . 5*V21 *VAR ( 1 ) 
PRINT 104,XK,CKSL3P,CKL4P, CKSM3P* CKM4P, ( VAR ( I ) ,1=2,9) 

WRITE (7,1 10 )XK,CKSL3P,CKL4P,CKSM3P,CKM4P 
I F f XS.LT.XSTOP)GO TO 9 $ ENDFILE 7 $ GO TO 1 
999 REWIND 7 
END 


PROGRAM PF1STRP 
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SUBROUTINE INICIDPDXR t OPDXI,XO,THW,BO,PRO,PID,TNL) 
************************************************************************* 

* SUBROUTINE INIC CALCULATES THE INITIAL VALUES AND DERIVATIVES FOR * 

* BEGINNING THE NUMERICAL INTEGRATION AT XO. * 

************************************************************************* 

C OM MON/BLOCK I /CUVAR (9) , HER (9 ) ,VAR{9) , PDB , RODB , UOB , VDR , XM2 S , XBO S , 

+ PBOS,CV$R,CPLGI,CPTHI,CTNL2,DVSDXR,DVSDXI,RVDB, PDV.UDV, VDV.RCNL 
+ , OR 0 , D 1 0 

VSR=CVSR $ VS I=CPLGI $ DVSDXI=CPTHI t DVSDXR=0. 

TNL=TAN( BO-THWI $ SNL=S I N < BC-THW ) $ CNL=COS ( BO-THW I $ RCNL=I./CNL 
BRO=(VSR+VDV*OIO)/VDB t PR 0= PDB *BRO-PDV*D 1 0 
BIO=( VSI-VDV*DRO)/VDB t P I 0= PDB*B 1 0 + PDV*DR0 
URG=UDR*8R0-UDV*DI0 i U I 0=UDB*B I O+UDV *DRO 
PRINT I ,VSR »VSI ,8R0,BI0,PR0, PI Q , UP.0 ,U 1 0 , DRO , D I 0 

1 FORMAT ( /* VSR=*G16 . 8 » * VSI=*G16.8,* BRO=*G Ifc . 8 , * BI0=*G16.8/ 

♦ * PRO=*G 16 . 8 » * PI0=*G16.8,* UR0=*G16.8,* UI 0=*G1 6 . 8/ . 

♦ * DR0=*G16.8,* DI0=*Gl6.8/> 

BNUM1=1.-XB0S*TNL*TNL t BNUM2= . 5*XB0S*TNL*TNL 
BNUM3=TNL* ( XM2S*PDV-UDV+ . 5 *VDV* X80S*TNL I 

B NU M4=TNL*( XM2S*PCB-UDB+V0V/ SNL* . 5*XB0S*TNL *< VDB + 2. *PDV/ SNL )) 
DENOM=VOB+XBOS*TNL*PDB 

D BD XP= ( BNUM 1 *DV SDXR + BNUM2* VS I + BNUMS^DRO+BNUMA *BI0)/DE NOM 
DBDXI=( BNUM1*OVSDXI-BNUM2*VSR+BNUM3*DIO-BNUM*»-*BRO)/DENOM 
DPDXR=PDB*DBDXR+TNL*(0VSDXR-.5*( VSI+VDV*DR0+BID*(VDB+2.*PDV/SNL) ) ) 
DPDXI=PDB*DBDXI*-TNL*(DVSDXI+.5*( VSR-VDV*D l 0+BR0*( VDB+ 2. *PDV/ SNL ) > I 
DUD XR=-0 PDXR + UI 0 $ DUOX I DPCX I-URO 
PRINT 2,DBDXR,DBDXI ,DPDXR,DPDXI ,DUDXR,DUDXI 

2 FORMAT!* DBDXR=*G16 .8,* D BDX I=*G 16 . 8, * DPDX R= *Gl 6 .8 , * DPDXI=*G16.8 

«■»/ * DUDXR=*G16. 8,* DUDX I=*G16 . 8/ ) 

VAR (2)=BR0*X0*0BDXR $ V AR ( 3 > =B I 0+X0*DBDX I $ VAR(1)=X0 

VAR (A) = PRO«-XO*OPDXR $ V AR ( 5) =P 1 0+X0*DPD XI 

VAR (6 )=URO+XO*DUOXR $ V AR ( 7 ) =U I 0+XQ*DUDX I 

VAR (8 1=DR0+XC*8R0*RCNL t V AR (9 > =D I 0+X0*B I 0*RC NL * RETURN 

END SUBROUTINE INIC 

SUBROUTINE CHSUB $ RETURN 

********* ************************* *************************************** 

* DUMMY SUBROUTINE CALLED BY INT2A. * 

************* ************************************************************ 

END SUBROUTINE CHSUB 
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SUBROUTINE DERSU8 

************** *********************************************************** 

* SUBROUTINE DERSUB CALCULATES THE FIRST ORDER DERIVATIVES FOR EACH * 

* STEP OF THE NUMERICAL INTERGR AT I ON . DERSUB IS CALLED BY INT2A. * 

************************** * * * **** * * * * * * * * * * * * * * * * * * $ * 4c * * * * * * * * * * * * * * * * * * * 

COMMON/BLOCK l /CUVAR (9) ,DER (9 » ,VAR(9) , PDB.RODB , UDB , VDB , XM2S. XBOS. 

+ RBOS»CVSR,CPLGI tCPTHItCTNL2,DVSDXR f DVSDXI , RVDB, PDV, UDV , VDV , RC NL 
+ , DR 0 , D I 0 

X= CUVAR ( 1 ) $ BR=CUVAR ( 2 ) t BI=CUVAR(3) $ PR=CUVAR(4) 

PI=CUVAR(5) $ UR=CU VAR { 6 ) $ UI=CUVAR(7) $ DR=CUVAR{8» $DI=CUVAR(9I 
VSR=CVSR t VSI=(CPLGI+X*CPTHI) $ RX=1./X 
DER f 2 )= ( (VDB*BR+CTNL2*( PR-PD B*BR ) + ( CTNL2*PDV- VDV ) *D I- VSR » *RX 
1 -DVSDXR+ ( ( VDB+VDV*RCNL J*BI+VSI+VDV*DR) )*RVDB 
DER (3 1= ( < VDB*BI+CTNL2*( P I- PD B*S 1 1- { CTNL2*PD V- VDV ) *DR- VS I ) *RX 
1 -DVSDXI- ( ( VDB+VDV*RCNL »*BR+VSR-VDV*DI) )*RVDB 
DER (4)= ( RX*(CTNL2*VSR+(CTNL2*VDV-XBOS*PDV)*DI+(PDB*XBOS-VDB* 

1 CTNL2>*BR-XBOS*PRI-XBOS*PDB*DER( 2) +< ( PDV4XM2 S-UDV) *DR-UI 

2 +XM2S*PI + ( PD B*XM2 S— UD8+XBOS *PDV*RCNL 1*81) ) *RBOS 

DER (5I=(RX*(CTNL2*VSI— ( CTNL2*VDV-XB0S*PDV ) *DR + ( PDB* XBOS- VDB* 

1 CTNL2J4BI— XBOS*PI )— XBOS*PCB*DER ( 3 ) +( ( PDV*XM2 S-UDV) *DI+UR 

2 -XM2S*PR— ( PDB*XM2S— UOB+XBOS*PDV*RCNL)*BR))*RBOS 
DER ( 6 ) = RX*( { PCB+UDB ) *BR-PR-UR- ( PDV+UDV )*DI )-DER (4) 

1 -( POB + UDB)*DER(2)+(UI + (UDB+(UOV + PDV)’*RCNL)*BI+UOV*DR) 

DER (7 )=RX*( ( PCB+UDB ) *BI-P I-UI+I PDV+UDV ) *DR I -DE P ( 5 » 

1 -( PDB+UDB) *DER (3)- (UR+ (UDB+ (UDV+PDV )*RCNL ) *BR-UDV*DI J 
DER ( 8 )= BR*RCNL $ DER (9)=BI *RCNL $ RETURN 
END SUBROUTINE DERSUB I 


SUBROUTINE SHKDF RV ( XM , G , THW , BO , ROO, UO , P V ,R OV, U V , VV, PB , ROB, UB , V B) 
************************************************************************* 

* CALCULATES SLOPE AND NORMAL VELOCITY DERIVATIVES OF P.RO.U, AND V FOR* 

* AN OBLIQUE SHOCK WAVE. DERIVATIVES ARE NORMALIZED 8Y| FLOW VARIABLES* 

* BEHIND THE SHOCK WAVE. * 

************************************************************************* 

SLO-SIN(BO-THW) $ C LO=C OS ( BO-THW J $ SB=SIN(BO) $ CB=COS(BO) 

SQMN=(XM*SB)**2 t FGI=4 . / ( G+ 1 . ) t RUO=l./UO 

PN= FGl*SB*RUO $ PV=PN/ROO $ PB=PV*CB*RUO 

RQV=FGI*ROO/( SQMN+SB) $ ROB=ROV*CB 

UN=;5*FG1*< SQMN+1. ) /SQMN $ UB=UN*RUO*SI N(THW)-PN*CLO 

UV=— UN^SLO $ VV=UN*CLO * VB=UN*RUO*COS( THW)-PN*SLO 

RETURN 

END SUBROUTINE SHKDERV : 
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SUBROUTINE WE DGE ( XM , G,THW, XM2, P20P1 , R20R 1 , V20V1 , BSTAO, I ER ) 
************************************************************************* 

* CALCULATES SHOCK WAVE ANGLE BET AO (RADIANS) FOR A WEDGE OF ANGLE * 

* THETA (RADIANS) FOR INPUT VALUES OF MACH NUMBER XM AND GAMMA G - * 

* FLOW CONDITIONS BEHIND THE SHOCK ARE CALCULATED FROM BETAO, XM»AND G* 
********************************** *************************************** 

R INDF=OI777COOCOOOOOOOOOD17 $ I F ( XM.GE . 1 . ) GO TO 1* IER=I $ GO TO 9 

1 IF(THW)8,7,2 

2 STSQ=SIN(THW)**2 $ SOM=XM*XM $ RSQM=1./SQM 

B=-RSQM*( SQM+2. )-G*STSQ $ D=( STSQ-1 . ) *RSOM**2 

C=< 2.*SQM+1 . ) *PSCM**2+( .25*( G* ( G+2 . >+I . ) +RSQM*( G-l. ) >*STSQ 

CALL WDETACH(XM,G»TDET) i IF(THW.LE.TDET)GO TO 3 t IER=3 $ GO TO 9 

3 X=-.333333333333333*B $ TCORR=I. $ DO A I=I,IOO 

CORR=-( D+X* (C+X*(B+X) t)/(C+X*(2.*B+3.*X>) $ XT=X t X=X+CORR 
I F ( I.NE. l.AND. (TCORR*CORR.LE.O. .OR.X.EO.XT) )G0 TO 5 

A TCORR=CORR t IER=A $ SMN=X*SQM $ GO TO 6 

5 I ER=C $ SMN=X*SQM 

6 P20 Pl = ( 2 • *G*S MN-G + 1 .)/(G + l.) $ R20R 1= ( G+l . ) *SMN/ ( 2 . «-SMN* ( G-I , ) ) 

V20V1=SQRT( 1 .—A . *( SMN-1. )* ( G*SMN+1 . ) / ( SM»*SQM* ( G*( G+2 . )+ I . ) ) ) 
XM2=XM*V20V1*SCRT( R20R1 /P20P1 I $ BETAO=ASI N ( SORT ( X) ) $ RETURN 

7 P20PI=R20R1=V20V1=1. $ XM2=XM $ BET AO=ASI N( 1 . / XM ) $ IER=0 $ RETURN 

8 I ER=2 

9 P2OPl=R20Rl=V20VI=BETA0=XM2=RINDF $ RETURN 

END SUBROUTINE WEDGE 

SUBROUTINE WDETACH( XM , G , DEL ) $ XMS=XM*XM $ GP=G+I. 
************************************************************************* 

* CALCULATES WEDGE ANGLE FOR SHOCK DETACHMENT DEL (RADIANS) FOR INPUT * 

* VALUES OF MACH NUMBER XM AND GAMMA G * 

************************************************************************* 

XMNS=.25*(GP*XMS-A.+SQRT(GP*(XMS*(GP*XMS+8.*(G-I. ))+I6 .) ) )/G 
S INDEL=SQRT ( ( XMS-XMNS* ( 2 .*XMS+1 ,-XMNS* ( XMS+2.-XMNS) ) ) / 

+ ( l.+XMNS*(G-l. + .25*XMS*GP*GP-G*XMNS) ) )/XM$ DEL = A SIN( SINDEL) ^RETURN 
END SUBROUTINE W DETACH 


I 
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SUBROUTINE INT2A( 1 1 ♦ N ,NT ,C I , SPEC, Cl MAX , IERR , VAR ,CUV AR , DER , ELE1 . 
1ELE2,ELT,£RRVAL,DERSU8,CHKSUB, I TEXT) 

***************** ******************************************************** 

* VERSION OF LRC CDC 6600 LIBRARY ROUTINE D2.4 * 

* USES FIXED INTERVAL, SINGLE PRECISION, 4TH-0RDER ADAMS-BASHFORTH * 

* PREDICTOR, 4TH-0RDER ADAMS- MOLL TON CORRECTOR, 4 TH- OR DER RUNGE- * 

* KUTTA STARTER. NT, C IMAX , IERR , ELT , AND ITEXT ARE DUMMY (UNUSED!. * 

********************************** *************************************** 

DIMENSION DER ( 21 ) * ELE 1 ( 20 ) , ELE2 ( 20 ) , EL T ( 20) ,ERR VAL( 20 ) , TEMP (20 I , 

1 DERl(20),DER2(20) ,DER3( 20) , S1VAR(21) ,VAR(2i) ,CUVAR(2I) 

A=0 .0 $ IF(II)1,1,2 

C INITIALIZATION SECTION 

1 IF(CI) 3,4,3 

C SAVE Cl 

4 PRINT 1000 * STOP 

1000 FORMATl 11H0CI=0 STOP) 

3 H=C I 

18 IERR=1 $ TO = VAR(l) $ M0DE=1 t 11= 1 S N1=N+1 $ DO 5 J=1,N1 
CUV AR(J ) = VAR(J) 

5 CONTINUE 

C EVALUATION SECTION HERE 

8 CALL DERSUB $ IF ( MODE .LE . 1 1 GO TO 7 * IF ( 1 1-3 ) 36 ,36 ,7 

36 CALL CHKSUB t IF ( 1 1 . EQ. 2» GO TO 1 

37 DO 38 J=1,N1 

38 VAR (JI = CUVAR(J) t I F ( 1 1-3 ) 6, 7 , 7 
7 RETURN 

6 I F( SPEC ) 9,7,9 ' 

9 OEL = VAR ( 1 ) -TO $ DELP=DEL* ( 1 .+ l.OE-11 ) 

I F ( ABS(OELP)— ABS(SPEC)) 2,10,10 

10 TO = VAR ( 1 ) $ GO TO 7 

2 11=1 $ I F ( MODE— 4 ) 11,12,12 

C RUNGE-KUTTA 

11 DO 20 J=2 , N 1 $ DER3(J— 1)=0ER2< J-l) $ DER2 ( J-l ) =DER1 ( J- 1 ) 

DER 1( J-l )=DER( J) $ ELE1 ( J-l )=DER( J ) $ CUVAR(J)=A 
DELT=0.4*ELEl(J-l)*H $ SI VAR ( J ) =VAR ( J ) $ CUVAR ( J )=S1VAR( J ) +DELT 

20 CONTINUE $ SI VAR ( 1 ) =V AR ( 1 ) $ CUVAR ( 1 )=S1VAR ( 1 )+0.4*H $ CALL DERSUB 


I 
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23 


24 

2 ? 


$ DO 24 J=2,N1 


I F ( 11-3)23,23,7 

CUVAR! 1>=S1VAR( 1 >+0.45573725421 879*H 

D EL T=f 0^296 97 7609 2477 5*ELE 1! J-l )+0. 15875964497104*ELE2( J-l ) >*H 
CUVAR(J)=SIVAR( JJ+DELT 

CONTINUE i CALL DERSUB $ I F( 1 1-3 ) 25 , 25, 7 
CUVAR ( 1 )=S1VAR( 1 l+H 1 $ 00 26 J=2,N1 $ T EMP ( J-l) = DE R! J ) 

OELT=(0 . 21810C38 822592*ELE 1( J- 1 >-3.05096514869293*ELE2( J-l ) 

1+3. 83286476046701*TEMP( J-l) I *H $ CUV AR I J ) =S IV AR { J ) +DELT 

26 CONTINUE $ CALL DERSUB $ I F ( 1 1-3 ) 27 , 27 , 7 

27 DH= H t CUVAR( 1)=VAR|(1 )+0H $ DO 28 J=2,N1 

OOUB= (0 . 1747602 82 26269 *ELE1 (J-l )-0 .5 514806628787 3*ELE2< J-l ) 
1+1. 2055 3 55° 93965 2*rEMP( J- l)+0. 1711847 8121 952* DER ( J) ) 

CUVAR (J) =VAR ( J ) +DH*|DOUB 
CONTINUE $ M00E=M0DE+1 $ GO TO 8 
AHAMS-MOULTON, ADAMS-BASHFORTH PREDICTOR 
CUV AR ( 1 ) =VAR ( 1 ) +H $! DH=H/24.0 S DO 13 J=2,N1 


28 


12 


13 

14 

15 


16 

19 


(55.0*DER( 


) =DER2 ( J) i DER2 ( J 
$ IF < 11-3)15, 15,7 


)=DER1 ( J) 


J)-59.0*0ER1 I J-l )+37.0*DER2< J- 1 ) -9 . 0*DER3 ( J-l ) ) 
VAR ( J ) + OH*|OOUB 

1 , N i DER3IJ 
t CALL DERSUB 
CTOR 

J = 2 » N 1 $ TEMfhCUVAR ( J) 

(9.0*DER( J)+19.*DERi( J-l )-5.0*DER2 ( J- 


D0U8 = 

CUVAR(J) 

CONTINUE $ DO 14 J= 
DERKJ ) = 0ER(J + 1) 
ADAMS— MOULTON CORR^ 
DO 16 
DOUB = 


CUVAR ( J )= VAR ( J)+DH4oOUR 


FRRVAU J-l) 
GO TO 8 
END 


MTEMP- 


• 1 ) +DER3 ( J- 1 ) ) 


CUVAR (J) ) /14. 21052631579847 


I 


SUBROUTINE INT2A 
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8.4.2 Program for Finite Difference Calculations 

The version of the program discribed here calculates either time 
asymtotic steady flows or flows with constant rates or accelerations. 
For an oscillation the flow field has to he converged for the initial 
rates and accelerations and the results used for the flow field for 
k = 1. Thus subroutine INIFLD would have to be modified slightly to 
perform an oscillation. 

8. 4.2.1 Input 

Each case consists of 11 cards of input data. Card 1 consists of 

80 columns of identification in an 8A10 format. Card 2 contains MSI, 

the number of (unstaggered) grid points between the shock and body; 

MFRES, the number of (unstaggered) grid points in the freestream on the 

initial line; LMAX, the number of grid points in the x-direction; and 

KMAX, the number of time steps in a 2014 format. Card 3 contains NFDTX, 

the number of CFL fractions; KMACH, the number of Mach numbers; and 

NTHtf, the number of wedge angles in a 2014 format. Card 5 contains 

FDTX, the CFL fraction array, in a 4E20.12 format. Card 6 contains of 

THETAD, the wedge angle array in degrees, in a 4E20.12 format. Card 7 

contains XO, the x for the initial line, and HXAMP, HYAMP, and THMAMP, 

the values of h^., hy, and 9 for oscillations, respectively in a 

4E20.12 format. Card 8 contains G, the value of 7 , in for ma t 4E20.12. 

Card 9 contains IOSC, which is 1 for this version, in a 2014 format. 

Card 10 contains HXDOT, HXDTDT, HYEOT, and HYETDT, the values of b^, 

• • • •• 

h x , hy, and h y respectively in format 4E20.12. Card 11 contains 
THMDOT and THMDTDT, the values of 0, and 0 respectively in a 4E20.12 
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format. When only part of the card is used, the leading fields are 
of course used. 

8. 4-. 2. 2 Output 

All case input parameters, effective wedge flow variables, the loop 
indexing for each line and k-step, and the flow variables on the initial 
line are printed. The version of the program listed here prints the 
flow variables for each body point for every 50 th time step and for the 


last four time steps. 
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8.4. 2. 3 Listing of Finite Difference Program . 


OV E RL AY ( L AX2 DU » 0 » 0 ) 

PROGRAM LAX2DU ( INPUT=1 ,OUTPUT=l,TAPE7=l, TAPE5=INPUT|^^^^^^^^ 

PROGRAM CALCULATES THE SUPERSONIC FLOW FIELD OVER A TWO-DIMENSIONAL * 
WEDGE AIRFOIL UNDERGOING SPECIFIED MOTION IN A PERFECT GAS. i? 

THE FINITE DIFFERENCE TECHNIQUE OF LAX IS USED. MAIN PROGRAM * 

PRIMARILY CONTROLS THE INOEXING AND LOOPING. SUBROUTINE INIFLD * 

READS THE CASE INPUT DATA, CALCULATES THE INDICES, AND CALCULATES * 
THE INITIAL FLOW FIELD. * 


* 

* 

* 

* 

* 

* 

♦#***#***j 


COMMON/ NLOOP/NFDTX,NMACHtNTHW 

COMMON/INDEX/NMEI 100) ,NMO( 100 ) ,NMTMOD(100 ) »MDIM, LDIM,LMAX»KMAX , 


+ 

1 


10 

20 

30 

50 

60 

70 

80 

90 

100 


200 

300 


MSI , MFRES 

NMACH=NTHW = NFOT X=1 

DO 300 MACH=1 , NMACH $ DO 300 NTHD=1,NTHW t DO 300 NFDT=1,NFDTX 
CALL INIFLO(MACH,NTHD,NFDT) 

DO ?00 K=2, KMAX $ KE VEN^MOD ( K ,2 > 

DO 100 L=2,LMAX $ LKBDY=M0D(L+K+1,2 > 

IFRES=MODINMTMOD(L) + LKBDYn,2) * MFL*NMO( D-IFRES 
IF ( KEVEN. EQ.O ) MFL=NME( L )— IFRES 

IFILKBDY.EQ.OIGO TO 50 $ IF( L.EQ.LMAXJGO TO 20 

CALL RFBDYPT (l,L+l»l,L,l,L-l*LKBDYJ $ DO 10 M=2 , MFL 
CALL LXFLDPT(M,L+1,M,L,M,L-1, M-l »L»K»LKBDY) $ GO TO 90 
CALL RFBDYPT(1,L,1,L-1,1,L,LKBDY> * DO 30 M=2 ,MFL 

CALL E XENDPT (M,L»M*L-1»M-1»L»K ,LKBDY) $ GO TO 90 

IFIL.EQ.LMAXIGO TO 70 * DO 60 M=1 , MFL 

CALL LXFLDPT(M,L+1,M+1,L,M,L-1,M,L,K,LKBDY) $ GO TO 90 
DO 80 M=1 * MFL 

CALL E XENDPT ( M+l ,L,M,L-i ,M,L,K,LKBDY) 

CALL MAXMPT (IFRES, MFL,L»K»LKBDY) 


CONTINUE 

CALL L SHI FT (KEVEN, K) 

IF(M0D(K,50).EQ.0.AND.K.LT. (KM AX- 3) (CALL FPRINT( KEVEN, Kl 
IF ( K.GE. ( KMAX— 3 ) I CALL FPRI NT ( KE VE N,K ) 

CONTINUE 

CONTINUE $ GO TO 1 

END PROGRAM LAX2DU 
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SUBROUTINE IN I FLO < MACH,NTHD, NFDT ) 

************************************* ****************************** ****** 

* THIS SUBROUTINE INITIALIZES THE ENTIRE FLOW FIELD FOR K=I AND * 

* CALCULATES AND STORES THE INDEXING PARAMETERS FOR THE STAGGERED * 

* GRID NETWORK AND FOR SUBSEQUENT TIME STEPS. * 

* INPUT OATA APE ALSO READ FOR MULTIPLE CASES. * 

************************************************************************* 

COMMON/NLOQP/NFOTX ,NMACH,NTHW 

COMMON/INDEX/NMEI 1001 ,NMO< 100 1 ,NMTM00( 100 » ,MDIM, LDIM,LMAX,KMAX, 

+ MSI »MFRES 

COMMON/FLOW PAR/P (25,100) »RO(25»100)*U(25»100),V(25»1 00) *RINDF » 

+ G.PF, SINT, COST 

CCMMON/MOT ION/ IGSC»X0, THMAMP , THMDOT , THMDTDT , SINTHM»COSTHM»HXAMP, 

+ HXDOT .HXOTDT « HYAMP , HYDOT,H YDTOT 
COMMON /GR I DPAR /DX , OT * OTX *GDT X,G12 

COMMON /FL INE / XM, THW, SI NOT, COSDT.P INI (25 ) , RO IN I ( 25 ) ,UIN I< 25 ) , 

+ VINK 25>,MFFI ,MSI2,MSI21 1 

DIMENSION IDENT(8I , DAT E( 2 > , AMACHI 4) ,FDTX (4 ) ,THETADI4) 

DATA MDIM,LDIM/25,100/ 

OATA FMT1,FMT2,FMT3,A10, 14 , E20/5H ( 1H1 > , 10H ( / X10A10/ ) ,6H(X914) , 

♦ 6H(8A10» ,6H(2014) ,9H( 4E20. 12 >/ 

DATA OEGRAD,TWOPI,RINOF/. 0174532925199433, 6. 2831 853071 796, 

+ OI77700000COQC0000013/ 

DO 1 J=1 , LOI M $ NME( J ) =NMO( J )=NMTMOD( J ) = 0 $ DO 1 I=l,MDIM 

1 P(I,J)=RO(I,J)=U(I,J)=V(I,J)=RINOF $ CALL DAYT IMCDATEI 
IF((MACH+NTHD+NF0TI.GT.3»G0 TO 20 

READ A10* I DENT * I F( EOF , 5 ) 999 , 2 

2 PRINT FMT1. 

PRINT FMT2,IDENT,CATE $ READ 1 4, MSI »MFRES,LMAX»KMAX 
READ I4,NF0TX, NMACH, NTHW $ READ E20, AMACH, FDTX , THETAD 
READ E 20, XO, HXAMP,HYAMP, THMAMPD, G 
READ 14, IOSC $ THMAMP=OEGRAO*THMAMPD 

IF(IOSC. EC. DREAD E20, HXDOT, HXDTDT,HYDOT ,HYDTDT, THMDOT .THMDTDT 
S I NTHM=0 • $ CO STFM=1 . 

20 IF ( IOSC.NE.OIGO TO 3 $ HXDTDT=HYDTDT=THMDT0T=THM=0. 

HXDOT=HXAMP $ HYDOT=HYAMP 

3 WRITE(7,AI0)IDENT ,DATE $ XM= AMACH( MACH ) $ DTXF=FDTX( NFDT ) 
THD=THETAO(NTHD) $ THW=DEGRAD*THD $ COST=COS( THW) $ SI NT=SI N( THW » 
VXE=COST+HXDOT * VYE=— SI NT+HYOOT 
EM=XM*SQRT(VXE*VXE+VYE*VYE ) t THE=ATAN(-VY E/VXEI 

CALL WEQGE(EM,G,THE,XMS,P20Pl,RS,US,BETA,IEWl $ I F, I I EW. NE.OI STOP 1 
PF=1/(G*XM*XM) $ PS=PF*P20P1$ TANBT=TAN( BETA-THE) S THEO=T HE/DE GRAD 
PRINT iOO,XM,G,THO,THE,DTXF,XMS, P20P1 ,RS,US,BETA,THED,EM, 

+ MOIM, LOIR, L MAX, MSI, MF RES 
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CFL1=1./!SQRT(2.>/XM+C0ST+SINT> $ CFL2=1 . / < SORT ( 2. )*US/XMS+US > 
CFL=AMINl(CFLltCFL2l $ OTXA=CFL*DTXF $ DX=XO*TANBT/I MSI-l.J 
OT=DTXA*DX $ IF! IOSC.NE.OIGQ TO A 

KMAX=2+THCPI/DT $ DT=T WO P I / ! KM AX-1 ) $ DTXA=CT/OX 

4 PRINT 104,XG,DX,DT,CFL,KMAX 

PRINT 105, IOSC,HXAMP f HXDOT,HXDTDT,HYAMP,HYDOT,HYDTDT , 

+ THMAMP.THMOOT.THMDTDT 

DTX=.5*DTXA $ GDTX=G*DTX $ G12=. 5*!G-1. ) t PRINT 1C1 
DC 6 L=1 * LM AX t NMS=MS 1 + 1 L-l » *TANBT 

NMT=NMS+MFRES $ NMTMOD ( LJ = MOD! NMT ,2 ) $ MGDL2=MOD! L, 2 ) 
NME<lJ=NMT/2+NHTMOD!U*MGD<L+l,2) t NMO! L) = NMT /2+NMTMOO ! L)*MODL2 
PRINT FMT3,L,NMO<t),NHE!L) ,NHT, NMTMOD (L) 

MSL=NMS/2+MQO(NMS»2) *MODL2 

DO 5 M=1 • MSL $ P!M,L)=PS $ RO!M,L)=RS $ UIH,L)-US 

5 V!M f LI=0. $ MFL=MSL+1 * MFF=NMO!L) $ DO 6 M=MFL.MFF 

P1M,L)=PF $ RO(M,L)=l.$ U! M,L)=VXE-I2*M-1-MODL2>*DX*THHDOT 

6 VIM,L)=VYE+IXa+DX*!L-l))*THMOOT 

MFFI=NMO!l) $ PRINT 102 $ 00 8 M=1,MFFI 

PINI !M»=P(M,1) $ RQINI !M l=RO !M» 1 ) % U IN I ( M » = U! M, 1 )$ VI NI ! M) =V i M, 1 ) 
8 PRINT 103»M»PINI!M) , ROINI (M) t UINI CM) » VI N I ( M ) 

MSI2=MSI/2 $ MSI21=MSI2+1 $ RETURN 

100 FORMAT!/* XM=*E21.13* G=*E21.13* THD,DEG.=*E21.13* THE=*E21.13 

1/* 0TXF=*E21 • 1 3* XMS=*E21.13* P2P01=*E21. 13/ 

2 * RS=*E21.13* US=*E21.13* BETA=*E21.13// 

3 * TH£D=*E21 .13* EM=*E21.13// 

4 * MOI M=* 1 4* LD IH=* 1 4* LMAX = *I4* MSI=*I4* MFRES=*I4//) 

101 FORMAT ( //3X*L NMO NME NMT NMTMOD* I j 

102 FORMAT !////3X*M I*9X*P INI *1TX*R0IN I *17X*U IN I*17X*VINI * ) 

103 F0RMATII4.4E21.13I 

104 FORMAT!/* XO=*E21.13* DX=*E21.13* DT=*E21.13* CFL=*E21.13* KMAX=* 
♦ I10//I 

105 FORMAT!* MOTION PARAMETERS AT K=1 , IGSC=*I4// 

1* HXAMP=*G21. 14* HXDOT =*G21 • 14* HXDTDT=*G21.14/ 

2* HYAMP=*G21. 14* HYD0T=*G21 . 14* HYDT0T=*G21.14/ 

3* THMAMP=*G21. 14* THMDOT =*G21»14* THMDT0T=*G21.14//) 

999 CONTINUE 

END SUBROUTINE! INIFLD 
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SUBROUTINE LXFLDPT(M1,L1 ,M2 , L2 ,M3 , L3,M4 , L4 , K .LKBDY) 

*************** ************************************************* ********* 

* SUBROUTINE CALCULATES FLOW VARIABLES AT A FIELD POINT GIVEN THE * 

* INDICES FOR THE FOUR BASE POINTS. ENTRY POINT BNDRYPT IS USED TO * 

* PERFORM S ARE OPERATIONS ON BOUNDARY POINTS GIVEN FLOW VARIABLES * 

* AT FOUR EFFECTIVE BASE POINTS. * 

************************************** **** *** ********* ****** ************* 

COMMON/FLCWPAR/P (2 5.100) ,R0( 25 ,100) , U< 2 5 .100 ) . V ( 25 » 100 ) « RINDF, 

+ G.PF, SINT, COST 

COMMCN/GRI CPAR/DX ,CT .CTX.GDTX.G12 

COMMON/MOT ION/ IOSC.XO.THMAMP , THMDOT , THMDTO T , S I NTHM ,COSTHM,HX AMP, 

+ HXDOT.HXDTDT.HYAMP.HYDOT.HYDTDT 

COMMON/GR IDVAR /PI ,P2,P3»P4»R01,R02» ROB »R04 ,U1»U2,U3,U4,V1,V2, V3, V4 
P1 = P (Ml , LI ) $ R01=R0( Ml, LI ) * U1=U(M1,L1) $ V1=V(M1,L1) 

P2=P ( M2, L2 ) $ R02=P0(M2,L2) $ U2=U(M2,L2) * V2=V(M2,L2) 

P3=P ( M3 »L3 ) $ P03=R0(M3,L3) * U3=U(M3,L3> $ V3=V(M3,L3) 

p4 = p ( M4 , L4 ) $ R04 = R0(M4,L4) $ U4=U(M4,L4» $ V4=V(M4,L4) 

ENTRY BNDPYPT 

RCU1=R01*U1 $ R0U2=R02*U2 * R0U3=R03*U3 $ R0U4=R04*U4 

ROV1=R01*V1 $ RCV2=R02*V2 * R0V3=R03*V3 $ R0V4=R04*V4 

R0US1=R0U1*U1 S ROUS2=ROU2*U2 t R0US3=R0U3*U3 $ ROUS4=ROU4*U4 

R0VS1=R0V1*V1 $ ROVS2=ROV2*V2 i ROVS3=ROV3*V3 $ R0VS4=RQV4*V4 

UROAV= . 25*( R0U1+RCU2+R0U 3+R0U4) $ ROAV=.25*(R01+R02+R03+R04) 

VROAV= .2 5* ( ROV1 + RQV2+ROV3+ROV4 ) $ DY=DX 

XLM=X0+DX*L3 * YLM=DY*( M2+M4+LKBDY-1 ) 

TRO=ROAV-CTX*( R0U1-R0U3+R0V2-R0V4 ) 

IF ( TRO* GT.O )G0 TO 1 $ IFD=1 $ GO TO 2 

1 RTR0=1 . /TRO 

TU=RTR0*{UP0AV-DTX*(P1-P3+R0US1-R0US3+R0U2*V2-R0U4*V4) 

+ -OT*<ROAV*(HYDTOT*SINTHM-HXDTOT*COSTHM)*THMDTDT*(ROAV*YLM 
+ +• 25*QY* (P02-P04) ) +THMCOT* ( 2 .*VR0AV-THMD0T*( ROAV*XLM 

+ *.25*DX*(R01-R03) ) ) ) ) 

TV=RTRO*( VROAV-OTX*(P2-P4+ROVS2-ROVS4+ROU1*V1-ROU3*V3) 

+ *OT*(ROAV*(HXOTOT*SINT HM+HYDTDT *COSTHM)+THMDTDT*(ROAV*XLM 
+ +.25*DX*(R01-R03) )*THMOCT*(2.*UROAV+THKDOT*(POAV*YLM 

+ +.25*0Y*(PC2-PC4)) ) ) ) 

TP=.25*(P1+P2*P3+P4)-GDTX*(P1*U1-P3*U3+P2*V2-P4*V4) 

+ +G1 2* ( ( . 25-DT X*U1 )* ( R0US1 +R0VS1 )>( . 25+DTX*U3 ) *( ROUS3+ROVS3 ) 

+ +( . 25-OT X*V2)*(R0US2+R0VS2)*(.25+0TX*V4)*( R0US4+R0VS4) 

+ -TRC* ( TU*TU+TV*TV) ) 

I F ( TP. GT. 0 ) GO TO 4 $ IFD=2 

2 PRINT 3, IFC ,K , LKBDY , Ml ,L1 , M2 ,L2, M3.L3 , M4.L4, PI ,P2, P3 »P4 , R01 ,R02, 

+ R03.R04 ,L1,U2,U3»U4,V1,V2,V3,V4»TP,TR0,TU,TV $ STOP 6 

3 FORMAT ( * IFD.K, LKBDY, ML 123 4= *1116/* PRUV1234T=*/(2X4E22.14/) ) 

4 P( M3 »L3 )=TP $ RO ( M3, L3 )=TRO $ U(M3,L3)=TU $ V(M3,L3)=TV $ RETURN 

END SUBROUTINE LXFLDPT 
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SUBROUTINE RFBOYPTI Ml ,L1 ,M2 , L2 , M3 , L3,K , LKBDY ) 
************************************************************************* 

* THIS SUBROUTINE ' TREATS THE BODY POINTS BY CALCULATING AN EFFECTIVE * 

* POINT INSIDE THE BODY USING THE REFLECTION APPROXIMATION WITH * 

* ALLOWANCE FOR DP/OY FROM UNSTEADY EFFECTS. * 

************************************************************************* 

COMMON/ FLOW PAR/P (25,100) ,R0 ( 25 , 10 0» ,U< 25, 100 ) , V( 25 , 100 I , RINDF , 

+ G,PF, SINT, COST 

COMMQN/GRI DPAR/CX*DT , CTX , GOT X « G1 2 

COMMON/MOTION/ IOSC , XO , THMAMP , THMDOT , THMOTOT ,SINTHM,CGSTHM»HXAMP, 

+ HXDOT »HXOTDT, HYAMP,HYDOT, HYDTDT 

COMMON /GR IDVAR/ PI , P2 , P3, PA , R01 , R02 ,R03 ,R04»U1*U2»U3*U4,V1*V2»V3*V4 
PI=PIM|,LI) $ R01=R0IM1,L1 ) $ UI=U(M1,L1) $ VI=V(MI,LI) 

P2=P(M2,L2) $ R02=RQ( M2»L2 ) t U2=U(M2,L2) t V2=V(M2,L2> 

P3=P(M3,L?) $ R03=RU( M3»L3 I * U3=U(M3,L3) $ V3=V(M3,L3) 

RG4=R02 $ U4=U2 $ V4=-V2 i M4=-M2 $ L4=-I $ DY =DX 

XLM=X0+(L2-1. )*0X 

P4'=P2~DY* (R01+RC3)*! HXOTOT*S INTHM+HYDT DT*COSTHM+XLM*THMDTDT 
+ + (U1+U3 )*THMOOT ) 

CALL BNDPYPT(M1,L1,M2,L2,M3,L3,M4,L4,K,LKBDY) $ RETURN 
ENO SUBROUTINE RFBOYPT 

SUBROUTINE MAX MPT (IFR £ S, MFL , L , K, LKBDY ) 
************************************************************************* 

* SUBROUTINE TREATS MAXIMUM Y-POI NT FOR A GIVEN L BY SETTING THE FLOW * 

* VARIABLES TO FREESTREAM VALUES OR TO INDEFINITES AS APPROPRIATE. * 
*************************************************** ********* ************* 

COMMON /FLOW PAR /P(25,100) ,R0! 25,100) »U( 25,100 ) *V(25*100) , RINDF, 

+ G,PF, SINT, COST 

COMMON /MOT I ON/ IOSC, XO , THMAMP , THMOQT ,THMDTDT,SINTHM,CCSTHM,HXAMP, 

♦ HXDOT, HXCTDT, HY AMP, HYDOT,HYOTDT 
COMMON/GRIOPAR /DX , DT , DTX , GOT X , Gl 2 

MF=MFLH * U=L-1 $ IFIIFRES.EQ.OJGO TO 1 $ XLM=X0*L1*0X 
YLM=(2*MFL-1-LKB0Y)*DX $ VX=HXD3T+C0ST $ VY=HYDOT-SI NT 
UBF=VX*COSTHM-VY*S I NTHM-YLM* THMDOT 
VBF=VX*SINTHM+VY*COSTHM+XLM*THMDOT 

PI MF ,L1 ) =PF $ ROI MF , L 1 ) = 1 • $ U(MF,L1)=UBF $ VIMF,L1)=VBF $ RETURN 
1 MIF = MF*IFRES $ P ( M I F , L 1) =R0( MI F , Li ) =U I MI F ,LI ) =V< MIF , H )=RINDF 
RETURN 
END 


SUBROUTINE MAXMPT 
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SUBROUTINE EXENOPTIM2.L2 ,M3, L3,M4,L4,K,LKBDY I 
****** ********** *************************************************** ****** 

* THIS SUBROUTINE EXTRAPOLATES THE AI TO DETERMINE THE FLOW VARIABLES 

* FOR CALCULATION OF THE LAST L COLUMN OF THE GRID. 
************************************************************************* 

COMMON/FLOWPAR/P I 25,1001 ,RO< 25,1001 ,UI 25 , 100 ) , VI 25 ,100 > , RINDF , 

+ G,PF, SINT, COST 

COMMON/GRI DPAR/CX,DT ,0TX,GDTX,G12 

C0MM0N/GRI0VAR/P1 ,P2,P3,P4,R01,R02,R03,R04,U1,U2,U3,U4,V1,V2,V3,V4 
P2=P (M2, L2) $ R02=RO( M2.L2 ) £ U2=U(M2,L2) $ V2=V(M2,L2) 

P3=P (M3, L3 ) $ R03=RQ( M3.L3) $ U3=UIM3,L3) $ V3=VCM3,L3) $ Ml=-1 

P4=P ( M4 , L4 ) £ RG4=R0(M4,L4 ) $ U4=U(M4,L4) $ V4=VIM4,L4) $ L1=L2+1 

R01 = R02«-R04-R03 £ IF I R01 .GT.O. )G0 TO 1 $ IEND=1 $ GO TO 2 

1 Ul= ( R02*U2+R04*U4 J -R03<‘U3 l/ROl £ V1=(R02*V2+R04*V4-R03*V3 )/R01 
PI =P2+P4— P3+G1 2 * ( 1102* ( V2*V2+U2* J2 ) +RQ4* ( V4*V4*U4*U4) 

♦ -R03*(V3*V3+U3*U:i)-Rai*<Vl*V2+Ul*Ul) I 
IFIP1.GT.0.) GO TO 4 $ IEND=2 

2 PRINT 3»IEND,K,LKIIDY»M2»L2»M3,L3»M4,L4,P1,P2»P3»P4»R01 ,R02,R03, 

♦ P04,U1»U2»U3,U4»V1»V2»V3,V4 £ STOP 10 

3 FORMAT ( * I END, K , LK8DY, ML 23 4= *9 1 6/* PRUV1234=*/ ( 2X4E22. 14/ )) 

4 CALL BNDRYPT(MI,L1,M2,L2,M3,L3,M4,L4,K,LKBDY) £ RETURN 

END SUBROUTINE EXENCPT 

SUBROUTINE L SH I FT ( KEVEN, K I 

************** ************ ***************************************** ****** 

* SUBROUTINE SHIFTS FLOW VARIABLES OVER ONE COLUMN IN THE TWO- * 

* DIMENSIONAL ARRAYSjANO CALLS INILINE TO DETERMINE NEW INITIAL LINE * 

* FOR L=l. * 

************************************************************************* 

COMMON/FLQWPAR/P ( 25,100) ,ROI 25,100) ,U< 25,1001 , V( 25 ,100 I , RINDF , 

+ G,PF, SINT, COST 

COMMON /I NDEX/NMEI 1001 , NMOI 1 00 1 ,NMTMOD ( IOC > , MDIM.LCIM.LMAX.KMAX, 

+ MS l , MFRES | 

00 1 M=1,MDIM $ P(M,LMAX)=RO(M,LMAX»=U(M,LMAX)=RINDF 

1 V(M,LMAXI=RINDF £ DO 2 L=2,LMAX £ LST=LMAX+2-L £ MFLT=NMO( LST ) +1 

IF(KEVEN.EQ.O) MFLT=NME <L ST )*i £ LL=LST-1 £ DO 2 M=l , MFLT 

P ( M» LST ) =P ( M » LL ) £ RO ( M, LST) =R0< M, LL ) £ U ( M , LST ) =U I M, LL ) 

2 VI M,LST)=V<M,LU £ CALL INILINE IKE VEN,K) £ RETURN 

END SUBROUTINE LSHIFT 


I 


* * 
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SUBROUTINE INILINEILKBDY »K) 

******************************** ******************* ********************** 

* THIS SUBROUTINE FILLS THE L = 1 LINE FOR EACH TIME STEP. FLOW * 

* VARIABLES ARE BASED ON THE QUASI-STATIC WEDGE ANGLE AND MACH * 

* NUMBER INCLUDING MOTION. CALLED BY L SHI FT ONLY. * 

*********************************** **************** *********** *********** 

COMMON /MOT ICN/ I OSC,XO » THMAMP ,THMDOT , THMDTDT ♦ S I NTHM»COSTHM,HX AMP, 

* HXDOT,HXDTDT,HYAMP,HYDOT,HYDTOT 

COMMON /F LC WPAR/P ( 25*1001 ,RO(25»iOO),U(25,100)*V(25»100) , R INDF , 

♦ G,PF, SINT, COST . 

COMMON/FLI NE/XM»THW*SI NOT , COSDT»PINI(25 ) .ROINl (251 ,UIN II 25 ) » 


♦ VINI(25I»MFFI*MSI2*MSI21 
COMMON /GR ICPAR/DXfDT ,DTX*GDTX»G12 

IF ( IOSC .EG .0 ) GO TO 2 * DO 1 M=1,MFFI 4 P ( M , 1 ) -PI NI (M ) 

RO ( M * 1 ) = RC INI ( M I $U ( M* 1 )=UIN 1 1 M ) 

1 V I M » 1 1 =V I N I < M 1 4 RETURN 

2 COSTIME=CCS(DT*(K-l. )) $ S INT IME= SINI DT* (K-l. )) 

THM=THMAMP*SI NT IME $ THMOOT=THMAMP*COST IME t THMOTDT=-THM 
S INTHM=SIN (THM ) $ COST HM=COS ( THM ) 

HYDOT=HYAMP*COST IME $ HYDTOT=-HYAMP*S INT IME 
HXOOT = HX AMP*CO ST IME $ HXDTDT=-HXAMP*SI NT I ME 
VXE=CQS (TFW-THM )-HYDOT*S INTHM+HXDOT*CQSTHM 
VYE=— SINITHW-THM) +HYOOT*COST HM*HXDOT*S I NTHM 
EM=XM*SQRT (VXE*VXE+VYE*VYE ) $ THE = AT ANI-VYE/VXE) 

CALL WE0GE(EM,G,THE,XMS,P20P1,RS,US,BETA,IEW) * IF( IEW.NE.OISTOP 3 
VX=HXOOT*COST $ V Y=HYDOT-S I NT $ UBFM=VX*COSTHM-VY*SINTHM 
VBF=VX*SINTHM+VY*COSTHM+XO*THMDOT 

DO 3 M=1 * MSI 2 $ P(M,1)=PF*P20P1 4 RC(M,1)=RS 4 U(M,1) = US 

3 VIM, 11=0. 4 DO 4 M=MS 121, MFFI 4 PIM,1)=PF 

R0(M,1)=1. 4 U(M,1)=URFM-DX*(2*M-1-LKBDY)*THMOOT 

A VI M, 1 ) =VBF 4 RETURN 

END SUBROUTINE INILINE 


SUBROUTINE FPR INT (KEVEN, K ) 

************************************************************************* 

* SUBROUTINE PRINTS INDEX AND FLOW VARIABLES FOR BODY POINTS AND * 

* ALSO WRITES THEM ON TAPE7 FOR USE 8Y PLOTTING PROGRAM. * 

***** ******************************************* ************************* 

COMMON/FLCWPAR/P I 25, 100>, RO( 25,100 ), UI 25,1 00 >,V! 25, 100) ,RINDF , 

♦ G,PF, SINT, COST 

COMMON/ 1 NDEX/NMEI100) , NMOI 100 ) .NMTMOOI IOC ) ,MDI M, LDIM , L MAX.KMAX , 

*■ MSI.MFRES 

PRINT 100, K 4 L1=2-KEVEN 4 DO 1 L=L1,LMAX,2 

IF (K. EQ.KMAX) WRITE (7,103 )L, PI 1,L) , ROl 1 , L ) , U < 1 , L ) , VII, L ) 

1 PRINT 101,L,P(1,L),R0U»L),U(1,L),V(1,L) 4 RETURN 

100 FORMAT!//* BODY POINTS FOR FINAL FLOW FIELD, L,P,FO,U,V, - K=*I6/I 

101 FORMAT (2X14, 4E 20.121 
103 FORMAT (I8«4E18.8) 

END SUBROUTINE FPR INT 
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SUBROUTINE WEDGE ( XM, G ,THW » XM2, P20PI * R20R1 , V20V1 » BET AO, IER) 
*********************************************************************** 

* CALCULATES SHOCK WAVE ANGLE BET AO (RADIANS I FOR A WEDGE OF ANGLE * 

* THETA (RADIANS) FOR INPUT VALUES OF MACH NUMBER XM AND GAMMA G - 

* FLOW CONDITIONS BEHIND THE SHOCK ARE CALCULATED FROM BETAO»XM,AND G* 


RINDF= 01777000 0000000000017 $ IF ( XM.GE. 1» ) GO TO 1$ IER=1 $ GO TO 9 

1 IFITHW 13*7,2 

2 STSQ=SIN(THW)**2 * SQM=XM*XM $ RSQM=1./SCM 

B=-RSQM*(SQM+2.)-G*STSQ $ D=< STSQ-1. )*RSQM**2 

C=(2.*SQM+1.)*RSQM**2+(,25*( G*( G+2 .)+l.)*RSCM*(G-l.) )*STSQ 

CALL WDETACH(XM*G,TDET ) $ IF (THW.LE.TOET )GO TO 3 $ IER=3 $ GO TO 9 

3 x=-. 333 3 3333333333 3* B 1 TCORR=l. $ DO A 1=1,100 

CORR=- (0 + X* (C+X* ( B+X ) ))/ (C+X*( 2. *B*3.*X ) ) t XT = X t X=X+CORR 
IF! I.NE.l.AND. ( TCORR*CORR. LE .0. . OR. X , EQ. XT I )GO TO 5 

A TCORR=CORR $ I ER=A * SMN=X*SQM $ GO TO 6 

5 IER=0 $ SMN=X*SQM 

6 P20P1=( 2.*G*SMN-G+1. ) /(G+l. ) 5 R20Rl=(G+l . ) *SMN/ ( 2.+SMN* (G-l . ) » 

V20Vl=SQRT ( l.-A.*( SMN-1. )*(G*SMN*1. ) / ( SMN*SQM* ( G* (G+2. )+l. )l ) 
XM2=XM*V2CV1*SGRT( R20R1/P2QP1 ) * BET AO=ASI N ( SQRT (X > ) S RETURN 

7 P20Pl=R20Pl=V20Vi=l. $ XM2=XM $ BETAO= ASI N ( 1. /XM ) $ IER=0 S RETURN 

8 I ER=2 

9 P20P1=R20R1=V2CV1=BETA0=XM2=RINDF $ RETURN 

END SUBROUTINE WEDGE 


SUBROUTINE WOET ACH ( XM, G, OEL ) $ XMS=XM*XM $ GP=G+1. 
************************************************************************* 

* CALCULATES WEDGE ANGLE FOR SHOCK DETACHMENT DEL (RADIANS) FOR INPUT * 

* VALUES OF MACH NUMBER XM AND GAMMA G * 

************************************************************************* 

XMNS = .25*(GP*XMS-A. + SQRT(GP*(XMS*(GP*XMS+8.*(G-1.) H-1S.) ))/G 

SINDEL=SQRT( ( XMS-XMNS* (2,*XMS+1.-XMNS*( XMS+2.-XMNS ) ) )/ 
+(1,+XMNS*(G-1. +,25*XMS*GP*GP-G*XMNS) ) )/XM$ OEL=ASIN( SI NDEL ) JRETURN 
end SUBROUTINE WOETACH 



